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Three  dynamical  features  associated  with  estuarine  lutoclines  are  examined  with 
special  reference  to  data  from  the  Jiaojiang  estuary  in  China.  These  features  include 
turbulence  damping  induced  by  high  suspended  sediment  concentration,  internal  wave 
behavior  at  the  lutocline,  and  the  response  of  the  lutocline  to  tidal  forcing.  It  is  shown  that 
turbulence  damping  is  governed  by  the  settling  flux,  cohesion,  interaction  between  floes  and 
sediment-induced  stratification.  Maximum  turbulence  damping  in  the  water  column  occurs 
at  the  lutocline,  a finding  which  supports  previous,  qualitative  observations  of  a similar 
nature.  Expressions  for  the  vertical  momentum  and  mass  diffusion  coefficients  incorporating 
these  effects  have  been  developed. 

Observations  from  the  Jiaojiang  are  examined  for  the  height,  angular  frequeney, 
celerity  and  length  of  internal  waves  at  the  lutocline.  Both  deep  water  high  and  shallow  water 
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low  frequency  waves  are  identified.  The  low  frequency,  at  0.09  rad  s , is  shown  to  be  close 


to  the  local  Brunt-Vaisala  frequency.  The  high  frequency  wave  at  1.33  rad  s is  possibly 

induced  by  interfacial  shear,  and  is  characterized  by  sharp  crests  and  flat  troughs.  The  height 
and  the  angular  frequency  of  both  wave  types  are  shown  to  decrease  with  increasing 
Richardson  number. 

Lutocline  variation  with  tide  in  the  Jiaojiang  is  examined  by  applications  of 
three-dimensional,  finite-difference  codes  developed  for  flow  and  sediment  transport.  It  is 
shown  that  lutocline  responses  reflect  the  cumulative  effects  of  sediment  settling  and 
entrainment,  turbulence  damping  and  tidal  asymmetry. 
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CHAPTER  1 
INTRODUCTION 

1.1  Problem  Statement 

A lutocline  is  a step  structure  in  the  vertical  profile  of  fine-grained  suspended 
sediment  concentration  (SSC),  and  is  a critically  important  pycnocline  in  addition  to  salinity 
and  temperature  gradients  governing  estuarine  dynamics  (Parker  and  Kirby,  1979).  From 
experiments  in  the  field  and  the  laboratory  in  recent  years,  the  lutocline’s  significance  in 
governing  the  vertical  mixing  of  suspended  sediment  and  sedimentation  patterns  has  been 
assessed  (e.g.,  Kirby,  1986;  Wolanski,  et  al.,  1989;  Mehta  and  Srinivas,  1993;  Winterwerp 
and  Kranenburg,  1997).  As  a result,  it  is  now  recognized  that  lutocline  formation  and 
strength  are  controlled  by  numerous  factors  including  tidal  currents,  interfacial  waves  formed 
at  the  lutocline,  turbulence  damping  in  the  fluid  mud  layer  beneath  the  lutocline,  floe  settling, 
suspended  sediment  deposition,  bottom  erosion,  and  consolidation  of  deposits. 
Unfortunately,  the  inter-linkage  among  these  unit  processes  has  yet  to  be  quantified  on 
theoretical  grounds  for  an  accurate  prediction  of  lutocline  dynamics.  There  is  a critical  need 
for  such  a quantification,  due  to  the  importance  of  predicting  fine-grained  sediment  transport 
in  estuarine  engineering  applications. 

Turbid  water  is  a characteristic  feature  of  estuarine  zones  with  high-load  fine 
sediments,  especially  around  turbidity-maxima,  sedimentary  fronts,  and  dredging  disposal 
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sites.  In  such  cases,  the  spatial  distribution  of  SSC  and  the  horizontal  transport  of  sediment 
are  significantly  associated  with  the  high  concentration  distribution  near  the  bottom,  i.e.,  the 
lutocline  and  fluid  mud.  Lutoclines  frequently  can  cause  significant  problems  due  to 
associated  sedimentation  in  engineering  projects,  especially  where  the  natural 
environment  is  dramatically  altered  by  dredging  operations  or  erection  of  structures. 

Prediction  of  lutocline  dynamics  is  presently  hindered  by  a lack  of  adequate 
understanding  of  lutocline-associated  processes,  namely  the  influence  of  SSC  on  turbulence 
damping,  hence  diffusion,  and  the  response  of  the  lutocline  interface  to  shear  flow  that  leads 
to  interfacial  waves  and  vertical  mixing. 

A high  degree  of  turbulence  damping  is  known  to  occur  in  the  fluid  mud  layer 
(Wolanski,  et,  al.,  1992;  Kranenburg  and  Winterwerp,  1997).  Hence,  vertical  diffusion  and 
mixing  at  the  lutocline  are  influenced  not  only  by  the  sediment-induced  buoyancy  force,  but 
also  by  turbulence  damping.  It  has  also  been  shown  that  the  characteristic,  upward- 
asymmetric,  sediment  entrainment  caused  by  the  instability  and  breaking  of  internal  waves 
at  the  lutocline  occurs  concurrently  with  high  turbulence  damping  (Jiang  and  Wolanski, 
1 998).  However,  there  remains  a lack  of  detailed  theoretical  analyses  and  direct  evidence  of 
turbulence  damping  because  of  the  difficulty  in  observing  turbulent  damping  both  in  the  field 
and  in  the  laboratory. 

Substantial  efforts  have  been  devoted  to  understanding  lutocline  response  to  tidal 
forcing  (e.g.,  Wolanski,  et  al.,  1988;  Smith  and  Kirby,  1989;  Costa  and  Mehta,  1990;  Dong, 
et  al.,  1997).  It  is  found  that  the  elevation  and  strength  (i.e.,  the  steepness  and  size  of  the  SSC 
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step  structure)  of  the  lutocline  vary  with  the  tidal  cycle  under  the  combined  effects  of  the 
strongly  current-dependent  vertical  and  horizontal  sediment  transport  processes.  However, 
a comprehensive  quantification  of  lutocline  response  to  tide  is  still  difficult  because  of 
limited  field  data  and  process-based  formulations  for  the  interactions  between  SSC  and  flow 
as  noted.  This  difficulty  in  quantification  and  the  consequent  need  to  examine  SSC-flow 
interaction  was  the  main  motivation  behind  the  present  study. 

1.2  Objectives 

In  accordance  with  the  above  discussion,  the  objectives  of  the  investigation  were  set 
as  follows: 

1.  To  examine  the  effects  of  high  SSC  on  the  turbulent  mixing  length  over  the 
estuarine  water  column,  and  develop  an  expression  for  the  vertical  diffusion  coefficient 
accounting  for  these  effects. 

2.  To  investigate  the  behavior  of  the  interfacial  waves  at  the  lutocline,  as  a basis  for 
an  improved  understanding  of  the  vertical  mixing  processes  at  the  interface. 

3.  To  develop  a numerical  code  for  estuarine  flow  and  fine  sediment  transport, 
incorporating  the  latest  unit  proeess  models  for  sediment  transport  including  erosion/ 
entrainment,  diffusion,  settling,  deposition  and  consolidation. 

4.  To  test  the  model  against  laboratory  data  and  analytical  solutions,  and  apply  the 
model  to  extensive  field  data  obtained  at  Jiaojiang  estuary  in  China,  as  a means  to  comment 
on  the  utility  of  the  model. 
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1.3  Tasks 

With  respect  to  the  above  objectives,  the  specific  tasks  established  to  conduct  the 
study  consist  of  the  following: 

1.  Review  of  Previous  Studies:  Previous  studies  on  sedimentary  process  in  fine 
sediment  dominated  estuarine  environments  are  reviewed,  thereby  providing  the  groundwork 
for  the  subsequent  tasks. 

2.  Theoretical  Analyses:  A analytical  model  for  vertical  diffusion  in  the  water  column 
is  developed  and  incorporated  in  the  numerical  model  for  the  flow  and  sediment  transport 
developed  to  examine  tide-induced  lutocline  dynamics  (see  below). 

3.  Flow  and  Sediment  Transport  Model:  A three-dimensional,  finite-difference 
scheme  numerical  model  code  is  developed.  This  code  consists  of  two  major  parts,  namely, 
a flow  model  (called  Coastal  and  Estuarine  Hydrodynamic  model  - University  of  Florida  or 
COHYD-UF)  and  a fine-grained  sediment  transport  model  (called  Coastal  and  Estuarine 
Cohesive  Sediment  model  - University  of  Florida  or  COSED-UF).  The  model  system  focuses 
on  fine  sediment  transport  processes  and  flow-sediment  coupling  due  to  sediment-induced 
stratification  and  turbulence  damping  in  the  fluid  mud  layer. 

4.  Model  Parameters:  All  parameters  defining  the  flow  and  sediment  transport 
processes  required  for  COHYD-UF  and  COSED-UF  are  determined  from  works  of  previous 
researchers,  model  calibrations  as  well  as  prototype  data  analysis. 

5.  Modeling  Tests:  For  validation,  COHYD-UF  and  COSED-UF  are  tested  against 
some  special  forms  of  the  governing  equations  having  exact  or  approximate  solutions,  as 

• well  as  laboratory  and  field  data.  The  sub-models  for  consolidation  of  deposits  and  interfacial 


5 


entrainment  are  tested  against  experimental  data  of  previous  researchers.  Further  tested  are 
the  ability  of  the  developed  model  in  predicting  local  scour  using  data  at  a pier  in  a river 
(Appendix  D),  and  sedimentation  behind  a piled  structure  using  flume  test  data  (Appendix 
E). 

6.  Field  Data  Analysis:  Through  analysis  of  field  data  from  the  Jiaojiang  estuary,  the 
following  tasks  are  carried  out:  (1)  Examination  of  the  behavior  of  internal  waves  at  the 
lutocline  in  conjunction  with  previous  analyses;  and  (2)  Examination  of  lutocline  and  fluid 
mud  response  to  tidal  forcing. 

1.4  Thesis  Outline 

This  study  is  presented  in  the  following  order.  Chapter  2 introduces  the  vertical 
structure  of  fine  sediment  suspension,  definitions  of  lutocline  and  fluid  mud  layers,  causes 
of  fluid  mud  and  lutocline  generation,  influence  of  lutocline  on  turbidity  transport,  and  a 
review  of  turbidity  dynamics  in  several  estuaries. 

Chapter  3 describes  model  formulations  including  the  governing  equations  and  the 
boundary  conditions  for  flow  and  sediment  transport,  reviews  of  previous  studies  on  flow- 
sediment  transport  processes,  and  the  respective  mathematical  formulas.  Lastly,  the  results 
of  model  tests  are  reported  for  hydrodynamics  and  sediment  transport. 

Analysis  of  field  data  and  the  results  for  the  behavior  of  internal  waves,  tides,  SSC, 
settling  velocity  and  erosion  rate  in  the  Jiaojiang  are  covered  in  Chapter  4.  Also  introduced 
in  this  chapter  are  previous  models  for  the  behavior  of  internal  wave  height,  angular 
frequency,  celerity  and  wave  length  at  haloclines  and  thermoclines.  These  models  are  used 
to  explain  the  behavior  of  internal  waves  detected  in  the  Jiaojiang. 
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Chapter  5 presents  an  analytieal  model  for  the  turbulent  mixing  length  and  the 
associated  diffusion  coefficient.  Data  from  the  Jiaojiang  are  used  to  examine  the  validity  of 
the  model. 

Chapter  6 summarizes  the  hydrodynamic  and  sedimentary  formulations  and 
parameters  for  numerical  model  application  to  the  Jiaojiang.  Simulations  of  lutocline  and 
fluid  mud  dynamics  in  the  Jiaojiang  and  their  comparisons  with  data  are  then  presented. 

Based  on  above  results,  major  conclusions  derived  from  this  study  are  summarized 
in  Chapter  7. 

Appendix  A provides  derivations  of  the  governing  equations  for  hydrodynamics  and 
sediment  transport  with  respect  to  o-transformation. 

Appendix  B describes  the  numerical  technique,  including  the  back-tracing  approach 
used  in  the  discretization  of  the  inertial  terms  in  hydrodynamic  and  sediment  transport 
equations,  and  the  pre-conditioned  conjugate  gradient  method  used  in  solving  the  finite 
differential  equation  for  the  water  surface  elevation. 

Appendix  C presents  the  analyzed  results  for  the  effects  of  temperature  on  cohesive 
sediment  settling  using  the  experimental  data  of  previous  researchers. 

Appendix  D presents  the  simulated  results  of  contraction  scour  at  the  Haldia  Pier,  in 
the  Hooghly  River  in  India  using  COHYD-UF,  and  a comparison  with  observations  in  the 
field. 

Appendix  E presents  the  simulated  results  of  sediment  deposition  pattern  behind  a 
piled  barrier  carried  out  in  a flume  test  and  comparison  with  experimental  observations. 


CHAPTER  2 

LUTOCLINES  IN  ESTUARIES 


2.1  Vertical  Structure  of  Fine  Sediment  Suspension  and  the  Lutocline 
A typical  instantaneous  sediment  concentration  profile  and  associated  velocity  profile 
is  schematically  shown  in  Figure  2.1,  as  might  be  observed  in  a macro-  or  meso-tidal 
estuarine  environment  with  high  fine-grained  sediment  loads  (Mehta,  1989;  1991a).  A 
noteworthy  characteristic  of  the  vertical  profile  of  the  suspended  sediment  concentration 
(SSC)  is  the  multiple  step  structures  identified  as  secondary  lutoclines  and  the  primary 
lutocline.  A secondary  lutocline  (Mehta  and  Li,  1997)  has  a relatively  low  SSC,  generally 
~1.0  kg  m'^,  and  is  induced  by  the  coupling  between  concentration  dependent  settling 

velocity  (i.e.,  settling  velocity  increasing  with  SSC  as  the  floes  become  stronger,  larger  and 
denser  in  the  flocculation  settling  mode;  see  Figure  2.2),  and  the  SSC  gradient  dependent 
diffusion  (i.e.,  the  upward  mass  diffusion  retarded  by  this  gradient  due  to  negative 
buoyancy).  The  primary  lutocline  is  a significant  pycnocline  which  occurs  near  the  bottom. 
It  is  characterized  by  a high  value  of  SSC,  generally  > -5-10  kg  m and  a significant 

vertical  gradient  of  SSC.  It  is  formed  due  to  a high  sediment  settling  flux  above  it  and  high 
turbulence  damping  as  well  as  hindered  settling;  see  Figure  2.2  below  (Ross  and  Mehta, 
1989;  Mehta  and  Li,  1997). 
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Figure  2.1.  Vertical  SSC  profile  classification  and  associated 
velocity  profile.  Also  shown  are  unit  transport  processes  which 
govern  concentration  profile  dynamics  (after  Mehta,  1989;  1991a). 
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As  shown  in  Figure  2.1,  the  primary  lutocline  tends  to  separate  the  water  column  into 
two  reasonably  well-defined  zones.  The  near-bottom,  non-Newtonian  flow  zone  is 
characterized  by  high  SSC,  generally  > ~5-10  kg  m'^,  hindered  settling,  turbulence 
damping,  turbidity  underflow,  a sediment-modified  velocity  profile,  etc.  On  the  other  hand, 
the  overlying  zone,  generally  having  SSC  < 1.0  kg  m and  exceeding  2-3  kg  m during 

extreme  energy  events,  is  marked  by  free  or  flocculation  settling  (Figure  2.2),  a higher  level 
of  turbulence  and  near-Newtonian  flow  properties. 

The  typical  SSC  profile  can  be  sub-divided  into  a layer  of  mobile  suspension,  the 
lutocline  layer,  mobile  fluid  mud,  stationary  fluid  mud,  a consolidating  mud  layer  and  the 
fully  consolidated  bed.  The  lutocline  layer  is  characterized  by  a thickness  6^,  between  the 

upper  elevation,  Z^,  and  lower  elevation,  Z^,  of  the  vertical  gradient  of  SSC  (Figure  2.1). 

The  characteristic  elevation  of  the  lutocline  above  the  bottom,  C , is  within  the  lutocline 

layer.  Approximately  at  elevation  Z^,  the  lutocline  layer  transitions  to  layers  consisting  of 

horizontally  mobile  and  horizontally  stationary  fluid  muds  (Figure  2.1).  Within  the  stationary 
fluid  mud  layer,  the  sediment  is  maintained  in  suspension  by  turbulent  diffusion.  Below  this 
is  the  consolidating  mud  layer. 

2.2  Causes  of  Fluid  Mud  and  Lutocline  Generation 
As  noted,  the  primary  lutocline  is  characterized  by  a steep  vertical  gradient  of  SSC 
and  occurs  within  a thin  transition  zone  (or  the  lutocline  layer)  from  the  upper  mobile 
suspension  to  the  lower  fluid  mud  (Ross,  1988).  Thus,  the  lutocline  occurs  concurrently  with 
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the  fluid  mud.  The  formation  of  fluid  mud  and  lutocline  are  dependent  on  complex  erosion, 
deposition  and  mixing-settling  processes,  high  turbulence  damping  due  to  settling  flux 
(Einstein  and  Chien,  1955),  cohesion,  interactions  between  floes  (Bagnold,  1954  and  1956), 
and  the  buoyancy  effect  due  to  sediment-induced  stratification  (Ross  and  Mehta,  1989).  A 
general  description  of  the  settling  behavior  as  well  as  causes  of  fluid  mud  and  lutocline 
generation  are  given  in  the  following  sections. 


Figure  2.2.  A representative  description  of  settling  velocity 
and  flux  variation  with  SSC  (after  Mehta  and  Li,  1997). 
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2.2.1  Settling  Velocity 

In  water  with  salinity  >0.1-0.5%o,  the  settling  velocity,  of  cohesive  sediment  is 

strongly  dependent  on  SSC  (=c),  and  can  be  divided  into  three  sub-ranges  in  terms  of 
concentration.  These  ranges  are  free  settling  (c<c, ),  flocculation  settling  (c,<c<C2)  and 

hindered-settling  (OCj)  with  the  maximum  SSC  for  free  settling  c, -0.1  -0.3  kg  m the 

maximum  SSC  for  flocculation  settling  Cj^S-lO  kg  m (Mehta  and  Li,  1997)  and  the 

maximum  SSC  for  hindered  settling  Cj =~80  kg  m (Odd  and  Rodger,  1986).  Free  settling 

occurs  at  low  SSC  when  the  settling  velocity  is  independent  of  concentration.  In  the 
flocculation  range,  the  settling  velocity  increases  with  concentration  due  to  the  formation  of 
stronger,  denser  and  larger  floes  with  increasing  concentration.  In  the  hindered  range,  the 
occurrence  of  an  aggregated  particulate  network  inhibits  the  upward  escape  of  the  interstitial 
water  present  within  the  deposit.  As  a result,  the  settling  velocity  decreases  with  increasing 
SSC.  Figure  2.2  gives  a representative  description  of  the  settling  velocity  and  the  settling  flux 
(=o)^c),  variation  with  SSC,  where  a,  b,  a and  P are  sediment-dependent  empirical 

coefficients  in  the  combined  relationship  between  the  settling  velocity  and  SSC  for  the 
flocculation  and  the  hindered  settling  regions,  is  the  free  settling  velocity,  is  the 

maximum  flocculation  settling  velocity,  and  is  the  maximum  settling  flux. 
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2.2.2  Formation  of  Fluid  Mud 

Fluid  mud  can  form  during  rapid  erosion  or  deposition.  During  erosion,  if  the  erosion 
rate  greatly  exceeds  the  net  upward  flux  due  to  sediment  diffusion  and  settling,  i.e.,  the  rate 
at  which  sediment  is  mixed  by  turbulence  into  the  upper  column  mobile  suspension  layer, 
the  SSC  near  the  bottom  continues  to  increase  and  tends  to  lead  to  the  formation  of  fluid  mud 
(Maa  and  Mehta,  1987).  During  deposition,  if  the  sediment  deposition  flux  exceeds  the  rate 
of  upward  transport  of  the  pore  fluid  (i.e.,  the  dewatering  rate  of  the  suspension),  a dense 
near-bed  suspension  is  formed,  grows  upward  and  consolidates  slowly  (Ross,  1988).  Once 
a fluid  mud  layer  forms,  it  is  stabilized  by  buoyancy  and  may  be  enhanced  by  high  turbulence 
damping  and  hindered  settling  as  well  as  further  erosion  or  deposition  (Maa  and  Mehta, 
1987;  Ross  and  Mehta,  1989). 

2.2.3  Formation  of  Lutocline 

Based  solely  on  the  concept  of  density  stratification,  the  lutocline  behaves,  at  least 
qualitatively,  like  a halocline  or  a thermocline.  Both  the  lutocline  and  other  pycnoclines  are 
characterized  by  the  coincidence  of  a vertical  maximum  in  the  concentration  gradient  and 
minimum  vertical  mixing  (Ross  and  Mehta,  1989).  The  local  minimum  in  mixing  due  to  the 
buoyancy  effect  is  associated  with  the  concentration  gradient.  However,  there  are  two 
significant  differences  between  a lutocline  and  other  pycnoclines  as  follows: 

(1)  Sediment  settling  occurs  with  respect  to  the  fluid.  This  has  the  effect  of  further 
enhancing  the  stability  of  the  lutocline  for  a given  set  of  flow  conditions  (Ross  and  Mehta, 
1989).  Based  on  the  non-linear  variation  of  the  settling  velocity  with  concentration  (Figure 
2.2),  at  the  maximum  settling  flux  the  gradient  ofdF ^dc  must  be  equal  to  zero,  i.e.. 


13 


dc  dc/dz 


(2.1) 


where  z is  the  vertical  coordinate.  Note  that  the  vertical  gradient  of  settling  flux  is  not  equal 
to  zero,  i.e.,  dF^dz*0,  due  to  the  high  gradient  of  SSC  at  the  lutocline.  Thus,  only  as 

dc/dz^°°,  Eq.  (2.1)  is  satisfied.  This  analysis  supports  the  observation  that  discontinuities 
in  the  concentration  profile  are  generated  at  elevations  where  the  maximum  settling  flux 

occurs. 

(2)  High  turbulence  damping  occurs  in  the  fluid  mud  layer  due  to  the  settling  flux, 
cohesion  and  interactions  between  floes  as  well  as  the  sediment-stratified  buoyancy  effect. 
Turbulence  damping  decreases  vertical  mixing  and  consequently  enhances  the  lutocline 
(discussed  in  details  in  Chapter  5). 

From  the  above  it  can  be  expected  that  the  lutocline  would  be  much  more  persistent 
in  the  macro-tidal  environments  than  the  halocline  and  the  thermocline  (Ross  and  Mehta, 
1989).  This  inference  has  been  confirmed  from  field  data  taken  in  the  Severn  estuary,  U.K. 
(Kirby,  1986),  with  a spring  tide  range  of  ~15  m,  and  in  the  Jiaojiang  (Dong,  et  al.,  1993, 
Guan,  et  al.,  1998  and  Jiang  and  Wolanski,  1998),  with  a spring  tide  range  of  6.3  m.  At  these 
sites,  lutoclines  have  been  observed  to  persist  through  most  of  the  tidal  cycle,  even  during 
spring  tides. 
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2.3  Influence  of  Lutocline  on  Turbidity  Transport 


In  a stratified  flow,  vertical  mixing  is  damped  due  to  the  buoyancy  force  because 
work  must  be  done  to  raise  heavy  fluid  parcels  and  lower  lighter  parcels.  As  noted,  in  the 
case  of  the  lutocline  the  mixing  process  is  further  complicated  by  sediment  settling,  cohesion 
(Hopfinger  and  Linden,  1982;  Ross  and  Mehta,  1989)  and  interactions  between  floes 
(Bagnold,  1954;  1956).  Thus,  additional  work  is  required  to  keep  the  sediment  in  suspension 
and  overcome  cohesion  and  interactions  between  floes. 

Without  explicitly  considering  settling,  cohesion  and  interactions  between  floes,  the 
mixing  process  over  the  lutocline  can  be  simply  demonstrated  with  reference  to  a sediment- 
induced  two-layer  flow  with  initial  densities  Pj  and  Pj,  and  identical  depths  //,  =H^=H/2\ 

this  system  representing  a fluid  mud  layer  beneath  clear  water  moving  parallel  to  the  jc- 
direction  with  velocities  t/,  and  U^,  respectively  (Figure  2.3).  After  some  time  complete 

mixing  is  assumed  to  take  place,  and  the  system  now  consists  of  single  layer  of  mean  density, 
p =(p,  +P2)/2  , flowing  with  mean  velocity,  U.  Because  initially  the  heavier  fluid  (of  density 

Pj)  lies  below  the  lighter  fluid  (of  density  Pj ),  the  center  of  gravity  of  the  system  is  below 

mid-depth  level,  whereas  in  the  final  state  it  is  exactly  at  mid-depth.  Thus,  the  center  of 
gravity  is  raised  in  the  mixing  process,  for  which  a potential  energy  PE  must  be  provided 
to  the  system  according  to 


0 


0 


(2.2) 
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where  g is  the  gravitational  acceleration.  The  source  of  potential  energy  increase,  PE,  must 
be  from  the  externally  imposed  kinetic  energy  so  long  as  the  initial  velocity  distribution  is 
uneven.  Conservation  of  linear  momentum  in  the  absence  of  the  effects  of  settling,  cohesion 
and  interactions  between  floes  implies  that  the  final,  uniform  velocity  is  the  mean  of  the 
initial  velocities  of  each  layer,  i.e.,  U=(U^  leading  to  a kinetic  energy  loss,  KE, 


given  by 


KE={ 


flp,uidz*fjp,ul 


dz 


H, 


H 


I Ip  U^dz~j-p(_U,-U,fH  (2.3) 


8 


Here,  the  Boussinesq  approximation,  p,=p2=p,  has  been  invoked  (Lesieur,  1997). 

Complete  vertical  mixing  is  possible  as  long  as  the  kinetic  energy  loss  exceeds  the 
potential  energy  gain,  i.e., 


(P2-Pi)g//^^ 

p{U,-U,f 


(2.4) 


Therefore,  the  initial  density  difference  must  be  sufficiently  weak  in  order  not  to  present  an 
insurmountable  gravitational  barrier,  or  alternatively,  the  initial  velocity-shear  induced  at  the 
interface  must  be  sufficiently  large  to  supply  the  necessary  amount  of  energy.  When  criterion 
(2.4)  is  not  met,  mixing  occurs  only  in  the  vicinity  of  the  initial  interface  and  cannot  extend 
over  the  entire  system.  This  situation  can  result  in  interfacial  instabilities  (Delisi  and  Corcos, 
1973).  This  phenomenon  will  be  further  discussed  in  Section  4.4.2. 
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Vertical  mixing  is  also  highly  dependent  on  the  nature  and  intensity  of  turbulence  in 
both  layers.  As  stated  in  Section  1.1,  an  upward-asymmetric  mixing  caused  by  the  instability 
and  breaking  of  internal  waves  at  the  lutocline  occurs  concurrently  with  high  turbulence 
damping  in  the  fluid  mud  layer. 


Figure  2.3.  Mixing  of  a two-layer  stratified 


Final  state  (after  complete  mixing) 
id,  with  fluid  mud  beneath  clear  water. 


2.4  Dynamics  of  Turbid  Estuaries 

Over  the  past  decades  there  has  been  an  increasing  number  of  studies  on  the  dynamics 
of  lutocline  and  fluid  mud  through  combined  physical  and  mathematical  approaches  and  in 
situ  measurements  of  sediment  movement.  This  has  been  made  possible  by  improved 
instrumentation  such  as  the  optical  backscatter  sensors,  e.g.,  turbidimeters  and  nephelometers 
(Wolanski,  et  al.,  1988),  and  acoustic  backscatter  sensors,  e.g.,  the  recently  developed 
Acoustic  Suspended  Sediment  Monitors  or  ASSM  (Shi,  1997).  Improvements  have  also 
occurred  in  designing  advanced  methods  of  recording  large  quantities  of  data  and  computer 
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analysis  techniques.  The  following  sections  hriefly  describe  the  works  of  previous 
researchers  in  some  of  the  turbid  estuaries  of  the  world. 

2.4.1  Amazon  Shelf 

The  Amazon  continental  shelf  in  Brazil  is  the  offshore  region  extending  from  the 
mouth  of  the  Amazon  River.  This  river  is  one  of  the  world’s  largest  with  the  highest  fresh 
water  outflow;  the  annual  mean  rate  of  about  210,000  m ^ s'*  is  1 8%  of  the  global  riverine 
discharge.  The  Amazon  also  has  one  of  the  largest  fluvial  sediment  discharge  of  about 
1.2x10^  tons  yr '* , or  10%  of  the  global  sediment  output  to  the  sea  (Kineke,  1993).  The 

shelf  region  can  be  considered  as  an  estuary  because  the  fresh  water  discharge  is  so  large  that 
low  salinities  (e.g.,  20%o)  occur  100  km  from  the  river  mouth  and  100  km  offshore. 

Measurements  of  SSC  and  sediment  transport  were  made  over  the  middle  and  inner 
Amazon  shelf  region  during  four  oceanographic  cruises  between  August  1 989  and  November 
1990  using  OBS  sensors  (range  0-36  kg  m CTD  probes,  current  meters  and  a data  logger 

mounted  on  a tripod  frame  (AmasSeds  Research  Group,  1990).  Silt-  and  clay-size  sediments 
were  found  to  dominate  the  bottom.  A major  characteristic  of  the  suspended  sediment  was 
the  vast  extent  of  high-concentration  near-bed  suspension  (or  fluid  mud),  which  accounted 
for  most  of  the  suspended  sediment  inventory  on  the  shelf.  Fluid  muds  generally  were 
observed  in  the  region  of  the  turbidity  maximum.  Fluid  mud  of  up  to  7 m in  thickness  on  the 
inner  and  middle  shelf  was  found  to  cover  an  extensive  region,  ranging  from  5,700  to  10,000 
km^,  in  water  depths  as  great  as  40  m.  Observed  SSC  near  the  surface  of  the  shelf  ranged 

from  0.01  to  >0.54  kg  m SSC  in  fluid  mud  layer  ranged  from  10  to  100  kg  m with  a 
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sharply  decrease  of  SSC  above,  thus  marking  a distinct  lutocline.  Much  of  the  inner  shelf 
mud  deposit  was  found  to  accumulate  sediment  at  a rate  exceeding  2cm  yr . 

2.4.2  Ariake  Bay 

Ariake  Bay  in  southern  Japan  is  dominated  by  semi-diurnal  tides  with  a spring  range 
of  3.86  m and  a neap  range  of  1 .54  m.  The  maximum  tidal  velocity  is  in  the  range  of  0.3-0. 5 
m s ■* . A great  amount  of  fine  sediment  is  transported  from  the  mountain  area  to  the  bay 

through  two  main  rivers,  which  lead  to  a wide  tidal  flat  near  Kumamoto  City.  The  bed 
material  offshore  consists  of  fine  clay  and  silt,  and  fine  sand  and  silt  occur  nearshore. 
Extensive  field  measurements  have  been  carried  out  at  Kumamoto  Port  in  order  to  obtain 
information  on  the  sedimentation  mechanism  in  the  port  area  (Tsuruya,  et  al.,  1990a). 

In  field  surveys,  an  echo-sounding  measuring  pole  (measuring  sedimentation 
volume),  electro-magnetic  current  meters  and  turbidity  meters  were  employed.  A noteworthy 
engineering  project  was  the  selection  of  a submerged  dike  for  reducing  deposition  within  the 
navigation  channel.  A multi-layered  numerical  model  was  used  to  take  into  account  the 
effects  of  the  dike  and  to  reproduce  the  vertical  distribution  of  suspended  mud  particles  in 
the  port  area  (Tsuruya  et  al.,  1990b).  It  was  found  that  there  was  a strong  correlation  between 
the  wave-induced  oscillatory  current  and  turbidity  concentration.  High  concentrations 
occurred  at  low  tide  when  wave  heights  were  large.  Maximum  concentrations  were  about  1 .0 
kg  m . SSC  greater  than  10.0  kg  m was  observed  near  the  bottom  with  a steep  vertical 

gradient  of  SSC  (or  lutocline).  It  was  also  found  that  net  bottom  erosion  took  place  when  the 
wave  height  was  large  and  tidal  level  low.  It  was  concluded  that  the  bottom  sediments  around 
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Kumamoto  Port  were  eroded  mainly  by  shear  stress  induced  by  wind  waves  (Tsuruya  et  al., 
1990a). 

2.4.3  Gironde  River 

The  Gironde  River  estuary  in  France  is  one  of  the  largest  estuaries  of  the  European 
Atlantic  coast.  It  has  an  annual  mean  inflow  rate  on  the  order  of  1 ,000  m ^ s',  and  the  tidal 

range  varies  from  2 to  6 m at  the  mouth.  The  range  increases  towards  the  upper  estuary  due 
to  the  convergence  effect  of  the  estuarine  shape,  and  tides  propagate  up  to  1 70  km  upstream 
during  low  river  inflow  periods.  Turbidity  maximum  is  an  important  feature  of  this  estuary. 
Numerous  studies  on  the  dynamics  of  this  turbidity  maximum  have  been  made  through  field 
investigations  (e.g.,  Allen,  et  al.,  1976;  Castaing  and  Allen,  1981),  as  well  as  depth-averaged 
and  three-dimensional  numerical  modeling  (Sottolichio,  et  al.,  1999). 

In  order  to  examine  the  resuspension  and  dispersion  of  fluid  mud  in  zone  of  the 
turbidity  maximum,  a tracer  study  was  carried  out  using  sediments  labeled  with  radioactive 
Scandium  46.  Five  kilos  of  labeled,  naturally  occurring,  fluid  mud  were  injected  into  a fluid 
mud  pool  in  a channel  during  a neap  tide  in  May,  1974.  It  was  shown  that  the  high  turbidity 
zone  was  characterized  by  SSC  ranging  between  1 to  10  kg  m (Allen  et  al.,  1976).  During 

slack  water  periods,  and  especially  at  neap  tides,  suspended  sediment  settling  was  enhanced, 
and  thick  patches  of  fluid  mud  appeared  on  the  channel  bottom  with  concentrations  up  to  300 
kg  m . The  estimated  total  sediment  mass  contained  in  the  turbidity  maximum-fluid  mud 

system  was  found  to  reach  up  to  5x10^  tons  (Castaing  and  Allen,  1981).  It  was  also 


determined  that  in  this  estuary,  tidal  pumping  is  the  major  mechanism  of  the  formation  of  the 
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turbidity  maximum  (Sottolichio,  et  al.,  1999).  Large-scale  lateral  transport  of  suspended 
sediment  occurs  within  the  turbidity  maximum,  resulting  in  a general  northward  drift  of  the 
sediment.  This  lateral  transport  occurs  by  advection,  while  the  rates  of  lateral  sediment 
diffusion  appear  to  be  low  (Allen,  et  al.,  1976). 

2.4.4  Hangzhou  Bav 

Hangzhou  Bay  on  the  coast  of  East  China  Sea  is  the  outer  region  of  the  Qiantang 
River  estuary  with  an  average  inflow  of  4.2x10*°  m^  yr  and  an  average  suspended 

sediment  load  of  7.9x  10°  ton  yr  '* . This  bay  is  a typical  funnel-shaped  estuary  with  a length 

of  about  100  km  and  an  average  depth  of  10  m.  Its  width  decreases  from  90  km  at  the  mouth 
to  20  km  at  its  western  end.  The  tidal  range  is  about  3 m at  the  mouth  and  increases  rapidly 
landward.  A tidal  bore  develops  about  10  km  further  upstream  of  the  bay  (Su,  et  al.,  1992). 
Sediment  in  the  bay  is  predominantly  fine  and  medium  silt  with  the  median  size  of  the 
suspended  load  ranging  from  1 0 to  1 3 pm  and  that  of  the  bed  about  1 6 pm  (Costa  and  Mehta, 
1990). 

Available  in  situ  data  are  from  the  following  three  surveys:  (1)  observations  in  the 
bay  carried  out  from  13'*  to  21'*  of  December,  1987  and  from  22'*  of  July  to  2'*  of  August, 

1988  using  Partech  7000-3RP  turbidimeter  and  ENDECO  174  current  meter,  aimed  at 
understanding  of  the  plume  front  and  its  role  in  suspended  sediment  transport  (Su,  et  al., 
1992);  (2)  Measurements  of  lutocline  and  flow-fine  sediment  hysteresis  near  the  south  shore 
were  conducted  between  14'*  and  16'*  of  May,  1988,  deploying  a pressure  gage,  two 


turbidimeters  (Partech  SDMI 6)  and  two  electromagnetic  current  meters  mounted  on  a tower 
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(Costa  and  Mehta  1990);  and  (3)  Fluid  mud  and  interfacial  waves  in  a proposed  navigational 
channel  observed  at  a neap  tide  (October  23,  1993)  using  a ship-borne  ASSM,  soon  after  a 
dredging  operation  (Shi,  1998).  Horizontal,  two-dimensional,  numerical  modeling  of  the 
depositional  patterns  in  Hangzhou  Bay  was  carried  out  by  Su  and  Xu  (1984).  It  was  found 
that  the  tidal  bore  there  traps  the  fine  fluvial  sediments.  A persistent,  year-around  NE-SW 
suspended  sediment  concentration  fi-ont  was  found  to  exist  inside  the  bay.  This  fi-ont  acts  to 
concentrate  suspended  sediment  and  transport  it  southwestward  into  the  bay  (Su,  et  al., 
1992).  ASSM  data  showed  that  high  concentration  suspensions  (SSC  > 10  kg  m “^) 

appeared  close  to  the  mud  bed,  occupying  30%  of  the  water  column.  The  ASSM  also 
detected  a relatively  well  defined  wave  train  with  a wave  period  of  about  144  s superimposed 
by  higher  frequency  oscillations  with  periods  of  a few  seconds  (Shi,  1998).  Due  to  the 
presence  of  the  lutocline,  reversals  of  flow-fine  sediment  hysteresis  was  observed  during  the 
transition  from  accelerating  to  decelerating  flows  (Costa  and  Mehta,  1990). 

2.4.5  James  River 

The  James  River  estuary  in  Virginia  discharges  into  the  south  end  of  the  Chesapeake 
Bay.  Extensive  surveys  of  the  vertical  profiles  of  SSC,  tidal  currents  and  sediment-water 
interface  were  conducted  using  an  optical  turbidimeter,  an  electromagnetic  current  meter  and 
a nuclear  transmission  density  probe,  all  mounted  on  a tripod  frame  (Nichols,  1984-1985). 
It  was  shown  that  the  flood  peak  at  6 cm  above  the  bed  reached  0.24  m s , whereas  the  ebb 

reached  0.30  m s Fluid  mud  accumulated  mainly  as  shallow  pools  and  blanket  deposits 


greater  than  0.2  m in  thickness  in  the  channel  depressions  in  the  middle  reach  of  the  estuary. 
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This  occurred  mainly  in  the  turbidity  maximum  zone,  a site  of  high  near-bed  concentration 
(0. 5-2.0  kg  m intensive  resuspension  and  rapid  sedimentation  (10-80  kg  m yr 

Measurements  of  mud  density  and  thickness  from  85  continuous  densitometer  profiles 
revealed  two  basic  types  of  profiles:  (1)  those  with  an  abrupt  increase  in  density  with  depth 
to  more  than  1,200  kg  m in  the  upper  1-2  cm,  i.e.,  mainly  settled  mud,  and  (2)  those  with 

a moderate  inerease  in  density  with  depth  in  the  density  range  -1,003-1,200  kg  m in  the 

upper  2-30  cm.  It  was  found  that  accumulation  of  fluid  mud  was  promoted  by  stratification 
of  the  interfacial  fluid  and  pore  water,  by  the  pseudoplastic  behavior  of  the  mud  with 
relatively  high  viscosity  at  low  shear  rates,  by  the  high  suspended  sediment  coneentrations, 
and  by  the  resultant  rapid-settling  flux  in  the  hindered  state  relative  to  the  consolidation  rate. 
2.4.6  Orinoco  River 

The  Orinoeo  River  delta  is  located  in  a coastal  plain  extending  from  the  Guyana 
shield  in  the  south  to  the  Venezuelian  corrdillera  in  the  north.  The  river  discharges  through 
the  delta  with  an  annual  mean  outflow  of  about  17,000  m ^ s and  total  suspended  load  on 

the  order  of  10*  tons  yr  “‘  (Eisma,  et  ah,  1978).  Bottom  deposits  sampled  in  the  delta 

(Eisma,  et  al.,  1978)  were  essentially  mixtures  of  silt  and  elay-size  materials,  with  generally 
low  sand  content  (less  than  1 0%). 

Turbidity  was  measured  using  the  radiometric  method  (Goldberg  and  Bruland,  1974). 
A layer  of  hyper-concentrated  suspended  material  (i.e.,  fluid  mud  layer)  with  concentrations 
reaching  500  kg  m was  located  at  the  bottom  of  the  navigation  channel  over  a distanee 
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of  up  to  50  km  and  a thickness  on  the  order  of  6 m.  This  layer  was  confined  to  the  navigation 
channel,  while  turbidity  exceeding  0.7kg  m was  not  detected  elsewhere  on  the  deltaic 

platform.  The  upper  layer  of  lower  turbidity  was  separated  from  the  fluid  mud  by  a sharp 
discontinuity  (lutocline),  which  caused  strong  reflections  on  the  echograms,  thus  resulting 
in  a so-called  “double-bottom”  pattern.  Within  the  fluid  mud  layer  there  was  a significant 
vertical  gradient  of  turbidity,  which  in  the  deepest  parts  of  the  channel  tended  to  reverse  close 
to  the  bottom  in  the  lowest  0.5  m of  the  water  column. 

2.4.7  Rhine-Meuse  River 

The  Rhine-Meuse  River  estuary  (or  Rotterdam  Waterway)  in  the  Netherlands  is  a 
partially  mixed  estuary,  through  which  the  Rhine  and  the  Meuse  flow  into  the  North  Sea  with 
an  annual  mean  discharge  of  1,500  m^  s . It  has  a deepened  navigation  channel  with  a 
depth  of  25  m.  In  field  investigations  of  the  movement  of  fluid  mud,  echo-sounders  and  a 
Partech  turbidity  meter  with  a range  of  0-12  kg  m were  employed  (Kirby  and  Parker, 
1977).  As  described  by  van  Leussen  and  van  Velzen  (1989),  fluid  mud  layers  with 
thicknesses  of  0.5-4.0  m and  SSC  < 15  kg  m usually  appear  in  this  area.  A distinct  lutocline 

was  observed  around  high  water  slack.  Fluid  mud  layers  play  a dominant  role  in  sedimentation 
in  this  estuary.  They  are  responsible  for  high  deposition  rates  in  short  times  under  rough 
weather  conditions.  At  times  sediment  deposition  of  about  2x10^  m ^ has  been  measured  in 


one  week  (Wiersma,  1984). 
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2.4.8  Severn  River 

The  Severn  River  estuary  in  the  UK  is  dominated  by  semi-diurnal  macro-tides  with 
a spring  range  of  13  m and  a neap  range  of  5 m.  Extensive  field  studies  in  this  estuary  were 
carried  out  in  the  early  1970's  (Kirby  and  Parker,  1982  and  Kirby,  1986).  In  these  studies, 
optical  turbidity  meters  were  used  to  obtain  continuous  horizontal  and  vertical  traverses  of 
SSC  in  comparatively  low  concentration  (0.1-20.0  kg  m'^)  areas,  and  gamma-ray 

densimeters  for  high-resolution  vertical  profiles  of  dense  stationary  suspensions.  The 
movement  of  fluid  mud  was  simulated  by  Odd  and  Cooper  (1989)  using  a horizontal,  two- 
dimensional,  numerical  model.  It  was  found  that  during  spring  tides  the  strong  tidal  currents, 
with  a flood  peak  on  the  order  of  2.3  m s'*  and  associated  turbulence,  were  sufficiently  high 

that  entrained  sediment  reached  the  water  surface  at  peak  ebb  and  peak  flood,  leading  to 
“well-mixed”,  homogeneous  vertical  SSC  profiles.  As  the  velocity  decreased  towards  slack 
water,  sediment  began  to  settle  whilst  the  velocity  was  still  high  (>1.5  m s'*).  A 

discontinuity  in  concentration  (i.e.,  a lutocline)  caused  by  settling  formed  in  the  water 
column,  which  then  subsided  towards  the  bed.  By  the  time  of  slack  water  the  layer  settled  to 
a level  only  2 or  3 meters  above  the  bed  to  leave  a residual  low  concentration  upper  zone  and 
a high  concentration  (~  >20.0  kg  m '^)  lower  layer.  At  slack  water  the  layer  stagnated  for 

a short  period  before  being  re-entrained  by  the  subsequent  reversal  of  the  tidal  current. 

2.4.9  South  Alligator  River 

The  South  Alligator  River  estuary  in  the  Northern  Territory  of  Australia  is  a macro- 
tidal  and  shallow  estuary  with  the  mean  water  depth  in  the  thalweg  on  the  order  of  only  a few 
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meters,  a spring  tidal  range  up  to  6 m and  maximum  spring  tidal  currents  of  up  to  2 m s . 

The  tidal  currents  exhibit  a strong  asymmetry  with  the  flood  currents  being  much  stronger 
than  ebb.  The  bed  sediment  is  a mixture  of  sand  and  mud.  The  bulk  of  the  suspended 
sediment  is  composed  of  particulates  1-4  pm  in  diameter.  Vertieal  profiles  of  SSC  were 
obtained  using  an  Analite  optical  fiber  nephelometer,  which  recorded  data  at  3 Hz  when 
lowered  at  a speed  of  about  0.2-0. 3 ms’*  through  the  water  column.  A 210-KHz  narrow- 

beam  Deso-10  acoustie  sounder  was  used  to  obtain  a visual  record  of  density  stratification 
induced  by  suspended  sediment  (Wolanski,  et  al.,  1988).  The  estuary  is  very  turbid,  with 
typical  SSC  values  on  order  1-6  kg  m . A maximum  SSC  of  10  kg  m was  measured. 

For  most  of  the  ebb  duration  a lutoeline  separated  a clear  upper  layer  from  an  extremely 
turbid  bottom  layer;  both  layers  being  of  comparable  thickness,  whereas  the  vertical  gradients 
in  SSC  were  small  at  flood  tides  (Wolanski,  et  al.,  1988).  A vertical,  one-dimensional, 
numerical  model  of  SSC  was  used  by  Wolanski,  et  al.  (1988).  Simulations  of  the  tidal 
evolution  of  SSC  were  shown  to  be  consistent  with  observations. 

2.4.10  Thames  River 

The  Thames  River  estuary  in  southern  England  is  characterized  by  relatively  uniform 
depths,  about  7.6  m at  mean  tide  level,  and  an  exponential  variation  of  the  cross  sectional 
area  and  charmel  width  along  its  length.  From  its  seaward  limit  to  the  tidal  limit  it  is  about 
100  km  long;  the  widths  at  these  limits  being  7000  m and  85  m,  respectively.  The  mean  tide 
range  at  the  mouth  is  4.3  m and  increases  up  to  5.6  m at  63  km.  With  one  major  exception, 
the  estuary  has  a hard  bed  made  up  of  gravel,  clay  and  chalk.  In  the  area  known  as  the  Mud 
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Reaches,  45-53  km  upstream  of  the  mouth,  there  are  extensive  deposits  of  silt,  and  a turbidity 
maximum  with  SSC  ranging  between  0. 1-5.5  kg  m ^^  (Owen  and  Odd,  1970).  Field 

investigations  (Inglis  and  Allen,  1957)  and  two-layer  numerical  modeling  of  Odd  and  Owen 
(1972)  found  that  the  null  point  is  usually  located  in  this  area,  and  the  turbidity  maximum 
and  sedimentation  there  are  due  mainly  to  silt  transported  in  the  lower  layers  jfrom  both  the 
upstream  and  downstream  directions.  A local  pocket  of  silt  also  occurs  at  about  24  km 
upstream  of  the  mouth. 

2.4.11  Yellow  River 

The  Yellow  River  estuary  in  northern  China  discharges  into  the  Bohai  Sea  with  an 
annual  mean  inflow  rate  of  about  1,550  m^  s"'  and  a sediment  load  of  about 

1.2x10^  tons  yr It  carries  the  largest  sediment  load  of  any  river  in  the  world  and  is 

dominated  by  loess  silt  and  fine  sand.  On  the  order  of  64%  of  the  sediment  is  deposited  on 
the  river  delta  and  mudflats,  while  the  remaining  is  transported  deeper  into  the  Bohai  Sea. 
As  a result,  a fan-shaped  delta  has  been  formed  since  1855,  which  covers  a distance  of  160 
km  along  shore  and  20-28  km  seaward  (Wang,  1988).  An  extensive  survey  over  the  active 
delta  front  was  carried  out  in  September-October  1987  and  July- August  1988,  aimed  at 
documenting  the  areal  extent,  vertical  thickness,  bulk  densities,  downslope  velocities  and 
velocity  gradients  of  the  sediment  underflows  and  their  relationship  to  tidal  and  storm 
forcing.  It  was  reported  that  cross-isobath  sediment  dispersal  into  the  shallow  Bohai  Sea  is 
dominated  by  the  formation  of  a hyperpycnal  plume  (SSC>2  kg  m ‘^)  and  gravity -driven 

underflows  (Wright,  1988).  The  strong  tidal  currents  often  had  speeds  of  over  1.5  m s 
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near  the  surface  and  0.8  ms'*  at  1 m above  the  bed.  The  observed  cross-isobath 
components  of  underflows  had  downslope  speeds  of  0.05-0.30  ms'*.  These  underflows 

descended  the  rapidly  prograding  delta  front  as  hyperpycnal  plumes  of  1-4  m thickness.  SSC 
in  the  lower  2 m of  the  water  column  in  the  shallow  parts  (<5  m)  normally  exceeded  1 
kg  m and  attained  maxima  of  over  10  kg  m '^ . SSC  near  the  surface  varied  from  <0.1 

kg  m to  > 1 .0  kg  m '^ . The  highest  turbidity  values  occurred  landward  of  the  front.  In  the 
deeper  part  (~10  m),  SSC  near  the  surface  was  consistently  less  than  0.1  kg  m and  about 
1 kg  m near  the  bottom.  Observations  during  a storm  and  immediately  after  the  storm 
revealed  near-bottom  layers  of  fluid  mud  with  a thickness  of  1-2  m and  average  SSC  at  -0.8 
m above  the  bed  of  about  252  kg  m '^ . Another  process  observed  at  the  top  of  the 

underflows  was  the  activity  of  low  frequency  (2.5-5.0x10'^  Hz)  internal  waves  of  relatively 

high  amplitude  (1.0-2. 5 m),  which  were  believed  having  contribution  to  plume  deceleration 
(Wright,  et  al.,  1988). 


CHAPTER  3 

FLOW  AND  SEDIMENT  TRANSPORT 
3.1  Introduction 

In  order  to  examine  lutocline  dynamics  in  estuaries  through  numerical  modeling,  it 
is  essential  to  carry  out  simulation  for  the  following:  (1)  ambient  flow  field,  (2)  flow- 
sediment  coupling  due  to  the  stratification  effect  and  turbulence  damping  in  the  fluid  mud 
layer,  (3)  erosion,  entrainment,  settling  and  deposition,  (4)  formation  of  fluid  mud,  (5) 
interfacial  mixing,  (6)  suspended  sediment  advection  and  diffusion,  (7)  consolidation  of 
bottom  sediment,  etc.  Accordingly,  a three-dimensional  numerical  model  code  is  developed. 
This  code  consists  of  two  parts,  namely,  a hydrodynamic  model  (called  Coastal  and  Estuarine 
Hydrodynamic  model  - University  of  Florida  or  COHYD-UF)  and  a fine-grained  sediment 
transport  model  (called  Coastal  and  Estuarine  Cohesive  Sediment  model  - University  of 
Florida  or  COSED-UF). 

This  chapter  begins  with  descriptions  of  the  governing  equations  of  hydrodynamics 
and  sediment  transport,  relevant  boundary  conditions,  fine  sediment  transport  processes, 
flow-sediment  coupling  and  finite-difference  schemes  for  solving  these  equations.  Results 
of  tests  related  to  model  validation  are  then  presented. 
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3.2  Hydrodynamics 

3.2.1  Goyeming  Equations 

In  the  deriyations  of  the  goyeming  hydrodynamic  equations,  the  following  treatments 
are  considered:  (1)  o-transform  is  introduced  (Figure  3.1)  by  transforming  the  temporal  and 
Cartesian  coordinate  system  (t,  x,y,  z)  to  a new  system  (/',  x',y',  o)  according  to  o^{z-(,)IH 

(Stansby  and  Lloyd,  1995),  (2)  the  flow  continuity  equation  takes  the  yertically  integrated 
form,  (3)  the  yertical  distribution  of  pressure  is  assumed  to  be  hydrostatic,  and  (4)  higher 
order  terms  related  to  diffusion  inyolying  o-coordinate  are  neglected.  Thus,  the  continuity 
and  momentum  equations  in  the  new  time  and  coordinate  system,  where  the  superscript 
(prime)  has  been  eliminated  for  conyenience,  respectiyely  are  (see  Appendix  A for 
deriyations): 


a 
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Figure  3.1.  Schematic  diagram  showing  o-transform. 
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Continuity  equation: 


dt 
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Momentum  equation  in  the  ^-direction: 
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Momentum  equation  in  the  y-direetion: 
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(3.3) 


where  D/Dt  denotes  the  total  derivative  with  respect  to  time,  u,  v and  o)  are  the  flow 
velocity  components  in  the  x,  y and  o directions,  respectively,  t is  time,  g is  the  gravitational 
acceleration,  C is  the  instantaneous  water  surface  elevation,  H is  the  total  water  depth, 
+h,h\s  the  imdisturbed  water  depth, /is  the  Coriolis  parameter,  A and  A are  the 

X y 
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horizontal  turbulent  momentum  diffusion  coefficients  in  the  x and>^  directions,  respectively, 

is  the  vertical  turbulent  momentum  diffusion  coefficient,  p is  the  fluid  density,  and  is  the 

Bingham  yield  strength  (Odd  and  Cooper,  1989).  The  vertical  velocity  in  the  o coordinate, 
(0,  is  determined  as  according  to 


dx  dy  J H 


dHu^dHv\^^ 
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(3.4) 


and  the  vertical  velocity  in  the  z coordinate,  w,  is  obtained  from  (see  Appendix  A) 
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3.2.2  Boundary  Conditions 

The  boundary  conditions  for  solving  the  above  equations  are  prescribed  as  follows: 
1 . Water  surface:  At  the  water  surface  (o=0),  the  x and  y components  of  the  wind- 

induced  stresses,  and  , are  respectively  specified  as 


X PA,du  V P^vdv 

where  the  resultant  vector,  f ^ , is  obtained  from 


(3.6) 


withC,=0.00l(l+0.07H)  (3.7) 

Here  p is  the  air  density,  C is  the  surface  drag  coefficient  and  W is  the  wind  velocity 


vector  at  a reference  elevation  (10  m above  the  water  surface  in  the  prototype  case). 


32 


2.  Bottom:  The  boundary  condition  at  the  bottom  (o=-l)  can  be  specified  by  either 
no-slip  or  a shear  stress  condition.  At  the  bottom,  the  water  particles  are  attached  to 
(stationary)  the  solid  surface,  hence  the  velocity  is  zero,  i.e.,  . Within  the  fluid,  velocity 


increases  rapidly  and  reaches  a value  of  (called  slip  velocity)  over  a small  distance 

(Woodruff,  1973).  Thus,  there  is  a very  steep  velocity  gradient  near  the  bottom.  If  one  were 
to  specify  the  no-slip  condition  right  at  the  bottom,  a high-resolution  numerical  grid  near  the 
bottom  would  become  essential.  In  turn  this  would  lead  to  overly  extensive  computation 
time.  In  addition,  viscous  effects  are  important  in  the  sublayer  near  the  bottom,  so  that  high 
Reynolds  number  turbulence  modeling  is  not  applicable  there.  Therefore,  the  shear  stress 
condition  is  most  commonly  used.  The  bottom  shear  stress  is  expressed  in  terms  of  the 
velocity  components  taken  from  the  modeled  layer  closest  to  the  bottom.  Accordingly,  the 

corresponding  stress  components,  and  x^,  can  be  related  to  the  velocity  gradient 
according  to 


_ P^v  du  ^_P'^vdv 


where  the  resultant  vector,  x^,  is  obtained  from 
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for  turbulent  flow.  Here  Zj  and  Fj , respectively,  are  the  elevation  above  the  bottom  and  the 


velocity  vector  of  the  modeled  layer  closest  to  the  bottom,  is  the  bottom  drag  coefficient. 
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K is  the  von  Karman  constant,  Zg  is  related  to  the  effective  roughness  of  the  bed  and  tj;  is  a 


stratification  function  which  takes  the  form  (Monin  and  Obukhov,  1953): 


'_z^ 

V 


= 1 +a, 


(3.10) 


In  Eq.  (3.10)  is  an  empirical  constant  in  the  range  of  4.7-5.S,  and  L^is  the  Monin- 


Obukhov  length  scale  defined  as 


Kgw'p'/pg  K^gz[(p^-p)/p]  dc/dz 


(3.11) 


where  D , is  the  bottom  friction  velocity  defined  as  w*=^|?j|/p,  c is  the  suspended  sediment 

concentration  (SSC),  Pg  is  the  density  of  water,  p^  is  the  sediment  granular  density,  w'  is 

the  turbulent  fluctuation  of  the  vertical  velocity  in  the  z coordinate  and  p'  is  the  turbulent 
fluctuation  of  density. 

Since  the  shear  stress  condition  is  employed  at  the  bottom  boundary,  the  flow  within 
low  turbulence  region  close  to  the  wall  can  only  be  described  by  a semi-empirical  wall 
function  which  bridges  the  viscous  sublayer  by  relating  the  values  at  the  first  numerical  grid 
point  placed  outside  the  viscous  sublayer  to  conditions  at  the  wall.  Launder  and  Spalding 
(1972)  have  proposed  a wall  function  which  is  described  by  a logarithmic  velocity  profile 
applicable  to  the  solid  wall  and  the  first  grid  point  adjacent  to  it.  The  standard  formulation 


of  the  wall  function  is 
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U 1 

— =-ln(y^€) 

K 


(3.12) 


where  y*=u^z^/\  is  a dimensionless  distance  or  frictional  Reynolds  number,  v is  the 


kinematic  viscosity  and  e is  a roughness  parameter  (e=9  for  a hydraulically  smooth  wall  and 
0.05  for  a hydraulically  rough  wall).  As  suggested  by  Rodi  (1980),  the  wall  function  should 
apply  to  a point  whose  y * value  is  in  the  range  of30<y  "^<100.  Eq.  (3.12)  is  then  sufficiently 
accurate  for  most  of  the  situations. 

4.  Shore  Boundaries:  At  the  shore  boundaries,  the  following  impermeable  condition 
is  prescribed: 


(3.13a) 

where  the  subscript  n,  denotes  the  direction  normal  to  the  solid  boundaries. 

Similar  to  the  bottom,  a no-slip  condition  at  the  shore  boundary  would  require  a high- 
resolution  grid  near  boundary.  Thus,  the  lateral  shear  stress  condition  is  specified  at  the  shore 
boundary  according  to 


dv, 


(3.13b) 


where  is  the  horizontal  turbulent  momentum  diffusion  coefficient  in  the  direction  normal 
to  the  solid  boundary,  is  the  coordinate  normal  to  the  solid  boundary,  is  the  tangential 
velocity  at  the  modeled  grid  closest  to  the  solid  boundary,  and  C,  is  the  lateral  friction 
coefficient. 


35 


5.  Open  Boundaries:  At  the  open  boundaries,  the  elevation  of  the  water  surface  is 
prescribed: 

Ce=Co+E^„cos(o)„r-g„)  (3.14) 

n = \ 


where  Cq  is  the  mean  water  level,  TV,  is  the  total  number  of  tidal  harmonic  components 

considered,  (o^  is  the  angular  frequency  of  the  «'*  harmonic  component,  and  and  are 

the  amplitude  and  phase  lag  of  the  harmonic  component,  respectively. 

3.3  Sediment  Transport 

3.3.1  Sediment  Conservation  Equation  in  the  Water  Column 

The  sediment  conservation  equation  in  the  and  o coordinate  system  can  be  stated 
as  (see  Appendix  A for  derivation  of  the  equation  in  the  o coordinate): 
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where  O)  is  the  sediment  settling  velocity,  and  AT,,  are  the  horizontal  turbulent  mass 

j X y 

diffusion  coefficients  in  the  x and  y directions,  respectively,  and  is  the  vertical  turbulent 


mass  diffusion  coefficient. 
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3.3.2  Boundary  Conditions 

The  boundary  conditions  for  solving  Eq.  (3.15)  are  prescribed  as  follows: 

1.  Water  Surface:  No  sediment  crosses  the  water  surface  (o=0),  i.e., 


1 c 1 1 


H da  Hda 


H da 


=0 


(3.16) 


2.  Bottom:  Sediment  vertical  flux  at  the  bed-fluid  interface  (o=-l)  is  equal  to  the  net 
flux  due  to  erosion  and  deposition,  i.e., 


1 1 1 


H da  Hda 


H da 


(3.17) 


where  w^is  the  flux  or  rate  of  bottom  sediment  erosion  (mass  per  unit  area  and  time)  and 
is  the  flux  or  rate  of  SSC  deposition. 

3.  Shore  Boundaries:  No  sediment  crosses  the  solid  shore  bovmdaries,  i.e.. 


dc 

dn. 


=0 


(3.18) 


4.  Open  Boundaries:  The  condition  at  the  open  boundaries  are  dependent  on  the  flow, 
i.e.,  inflow  or  outflow.  If  fluid  flows  into  the  modeled  region,  the  SSC  of  the  water  from  the 
outside  region  is  prescribed,  i.e., 

c,=cfx^,a,t)  (3.19) 

where  is  the  prescribed  SSC  at  the  open  boundary  and  c.  is  the  SSC  from  outside  the 

modeled  domain.  In  the  case  of  outflow,  the  diffusion  terms  in  Eq.  (3.15)  are  neglected  for 
numerical  purposes,  so  that  the  mass  conservative  condition  becomes: 
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(3.20) 


dt  dx  dy  da  H do 


3.3.3  Fine  Sediment  Transport  Processes 

To  use  Eq.  (3.15),  the  fine  sediment  transport  processes  must  be  expressed  in  terms 
of  mathematical  formulas.  The  major  relevant  processes  are  as  follows. 

1 . Consolidation:  In  cohesive  sediment  transport  models  there  is  a need  to  take  into 
account  consolidation  of  newly  deposited  sediment,  since  the  shear  strength  and,  thus,  the 
erosion  resistance  of  the  deposit  increases  with  consolidation  time  (Verreet  and  Berlamont, 
1 989).  The  erodibility  of  a mud  bed  is  determined  by  the  strength  of  its  particle-to-particle 
structure.  This  in  turn  is  influenced  by  the  state  or  degree  of  consolidation  of  the  bed,  which 
is  time-dependent.  Thus,  in  a consolidating  bed  the  resistance  to  erosion  is  a time-dependent 
function  of  consolidation  duration.  Consequently,  in  a model  of  fine  sediment  transport 
behavior,  predicting  bed  consolidation  occupies  a crucial  role. 

Consolidation  of  the  deposit  is  caused  by  the  self-weight  of  sediment  particles,  and 
is  a process  of  expulsion  of  pore  water  from  the  bed  (Parker  and  Lee,  1979).  The  settling 
cohesive  sediment  begins  to  behave  as  a soil  and  consolidate  when  the  resulting  stationary 
suspension  develops  particle-to-particle  contacts  and  an  effective  stress  appears  (Hayter, 
1983).  At  this  stage,  the  total  stress  within  the  soil  matrix  is  expressed  as 


where  is  the  total  (normal)  stress,  r'is  the  effective  (normal)  stress  and  is  the  pore 


(3.21) 


water  pressure. 
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The  aim  of  a consolidation  model  is  to  predict  the  vertical  distribution  of  the  time- 
dependent  dry  density  or  concentration,  c,  which  in  turn  is  related  to  the  soil  shear  strength. 
There  are  two  typical  theoretical  models  describing  the  consolidation  process,  those  based 
on  force  balance  (Gibson,  et  al.,  1967)  and  those  based  on  solid  mass  conservation  (Kynch, 
1952).  Here  the  theoretical  framework  first  developed  by  Kynch  (1952)  is  applied.  It 
describes  consolidation  by  a vertical  transport  equation.  The  consolidation  rate,  instead  of 
permeability,  is  used  as  the  model  parameter,  which  is  a function  of  dry  density  (Toorman 
and  Berlamont,  1993).  In  the  present  work,  the  difference  from  the  pure  deposition  model 
of  Kynch  (1952)  is  that  the  deposition  flux  minus  erosion  flux  at  the  bed-fluid  interface,  i.e., 
q=m^-m^,  is  incorporated  in  the  model,  and  a normalized  coordinate  o'  {=z'/H')  is  used. 


Thus,  the  relevant  transport  equation  is: 

dH'c 

= +q 

dt  do' 


(3.22) 


where  z'  is  the  vertical  coordinate  of  the  consolidating  layer  originating  from  the  bottom  and 
positive  upward,  i/'is  the  thickness  of  the  consolidating  layer  and  is  the  rate  of 


consolidation.  As  for  this  rate,  experimental  evidence  of  Toorman  and  Berlamont  (1993) 
suggests  that  two  distinct  modes  of  consolidation  can  be  recognized  from  the  plot  of 
consolidation  rate  against  concentration,  i.e.,  loose  soil  consolidation  and  compacted  soil 
consolidation.  For  these  two  modes,  they  developed  a combined  relationship  expressed  as 


( \ 

( \ 

ntt 

( \ 

w,iexp 

c 

F,+(o  , 

t sc2 

1-^ 

(l-F,),  F,=exp 

- 

c 

(3.23) 
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where  is  a characteristic  mode  transition  (from  loose  soil  to  compacted  soil)  function  with 


« >10,  m^  is  a sediment-dependent  constant  is  the  transition  concentration,  and  0)^^2 


are  the  rates  of  consolidation  for  the  first  and  the  second  modes,  respectively,  c^j  is  the 


concentration  corresponding  to  the  maximum  settling  flux  and  the  saturation 

concentration,  i.e.,  the  maximum  compaction  concentration. 

2.  Fully  Consolidated  Bed:  The  deposited  sediment  is  said  to  be  fully  consolidated 
when  c^c^2-  The  vertical  profile  of  the  dry  density  or  concentration  of  a fully  consolidated 


bed  can  be  expressed  as  (Mehta,  et  al.,  1982): 


c=Ca. 


f U )P2 
h.-LZ 


b ) 


(3.24) 


where  h,  is  the  bed  thickness,  C is  the  bed  average  concentration,  az  is  the  incremental 


depth  downward  from  the  bed  surface,  and  and  P2  are  bed-dependent  coefficients. 

3.  Shear  Strength:  The  bed  shear  strength  with  respect  to  erosion  is  the  primary 
measure  of  bed  scour.  As  stated  above,  the  shear  strength  increases  with  the  consolidation 
time,  or  the  dry  density.  Hence,  usually,  the  shear  strength  is  related  to  the  bottom  sediment 
concentration.  In  accordance  with  the  description  of  Mehta  (1991b),  the  bed  shear  strength, 
T , can  be  considered  to  have  the  form 

S ’ 


(3.25) 
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where  is  the  shear  strength  of  newly  deposited  bottom  sediment,  and  Pj  are  sediment- 

dependent  coefficients,  (j)  is  the  solids  weight  fraction,  (|)=c/p^,  and  (j)^is  the  critical  solids 
weight  fraction  below  which  mud  has  a fluid-like  consistency.  James,  et  al.  (1 988)  show  that  (j)^ 
is  typically  on  the  order  of  0.03  to  0.05. 

4.  Interfacial  Entrainment:  Since  diffusion-induced  mixing  over  the  lutocline  is 
damped  due  to  strong  stratification  and  turbulence  damping  within  fluid  mud,  internal  wave 
breaking  becomes  a major  mechanism  contributing  to  vertical  entrainment  over  the  lutocline 
(Scarlatos  and  Mehta,  1993).  There  are  two  primary  modes  of  instability  of  the  interface 
depending  on  the  relative  thicknesses  and  positions  of  the  current  shear  layer  of  thickness  6^ 

(Figure  2.1)  and  the  density  interfacial  layer  of  thickness  6^  (Mehta  and  Srinivas,  1993).  In 

case  the  mid-axes  of  density  and  velocity  gradients  coincide  and  6^  is  approximately  equal 

to  6^,  the  primary  mode  of  instability  is  of  the  Kelvin-Helmholtz  type,  and  is  characterized 

by  a roll-up  and  pairing  of  the  interfacial  vortices  (Delisi  and  Corcos,  1973).  When  6^  is 

smaller  than  6^  due  to  stratification,  Holmboe  type  of  instability  results,  as  recognized  by 

sharp-crested  interfacial  cusps  which  protrude  alternatively  into  both  fluids  (Browand  and 
Wang,  1972).  Results  from  laboratory  and  field  observations  further  show  that  turbulence 
damping  in  the  fluid  mud  layer  introduces  an  upward-asymmetric  mixing  over  the  lutocline, 
i.e.,  there  is  a net  upward  flux  of  mass  over  the  lutocline  (Wolanski,  et  al.,  1989;  1992; 
Mehta  and  Srinivas,  1993;  Kranenburg  and  Winterwerp,  1997;  Winterwerp  and  Kranenburg, 


1997;  Jiang  and  Wolanski,  1998). 
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Mehta  and  Srinivas  (1993)  established  a semi-empirical  formula  for  the  rate  of 
interfacial  entrainment,  which  accounts  for  the  cumulative  effects  of  settling,  cohesion 

and  viscosity  difference  (between  fluid  mud  and  water)  on  mud  entrainment 


^en  with  Ri^  - - 


(3.26) 


where  is  defined  as  dh^dt,  is  the  fluid  mud  layer  thickness,  and  are 
sediment-dependent  constants,  Ri^  is  the  global  Richardson  number,  C,  is  the  depth-mean 
sediment  concentration  in  the  upper,  mixed,  layer  of  height  is  the  corresponding 

value  for  the  lower  layer  (being  entrained),  (/,  and  (/,  are  the  respective  depth-mean  flow 


velocities,  and  is  a coefficient  dependent  on  the  granular  density  as 


m 


P.-Po 


(3.27) 


5.  Deposition:  From  flume  experiments  Krone  (1962)  concluded  that  the  rate  of 
deposition  is  equal  to  the  product  of  the  near-bed  settling  velocity,  SSC  and  the  probability 
that  a settling  floe  becomes  attached  to  the  bed: 


1--^ 


''dj 


(3.28) 


where  w . and  c,  are  the  near-bed  settling  velocity  and  SSC,  respectively,  and  t . is  the 


critical  shear  stress  for  deposition. 
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6.  Bottom  Erosion:  Erosion  of  cohesive  sediment,  which  is  dependent  on  the 
composition  and  structure  of  bottom  material  that  characterizes  bottom  resistance,  and  on  the 
nature  of  the  eroding  force,  can  occur  in  two  typical  ways  in  estuaries  (Mehta,  1991b).  The 
first  mode  is  floc-by-floc  surface  erosion  in  which  the  floe  at  the  bed-water  interface,  initially 
attached  to  their  neighbors  by  inter-particle  electro-chemical  bonds,  breaks  up  and  is 
entrained  as  a result  of  hydrodynamic  lift  and  drag.  The  second  mode  is  referred  to  as  mass 
erosion,  wherein  the  bed  fails  at  a deeply  embedded  plane  such  that  all  the  material  above 
that  plane  is  rapidly  brought  into  suspension. 

Surface  erosion  under  current-induced  bottom  stress  has  been  studied  extensively 
(Parchure  and  Mehta,  1985).  This  process  was  subsequently  examined  ftirther  by  Lee  and 
Mehta  (1994)  and  Mehta  and  Parchure  (1999).  From  these  studies,  the  effects  of  shear 
strength  and  temperature  on  the  erosion  rate  are  incorporated  as  follows: 

[0exp(A-A/0)]  (3.29) 


in  which  is  the  maximum  erosion  rate  constant  at  tj=2T^,  %,  X,  A and  A are  sediment- 

dependent  coefficients  and  0 is  the  absolute  temperature. 

7.  Settling  Velocity:  As  stated  in  Section  2.2.1,  in  water  with  the  salinity 
.s>0.1  -0.5%,  the  settling  velocity  of  cohesive  sediment  is  strongly  dependent  on  SSC  and 

can  be  divided  into  three  sub-ranges  in  terms  of  concentration.  These  sub-ranges  include  free 
settling,  flocculation  settling  and  hindered-settling  (Figure  2.2).  Hwang  (1989)  developed  a 
combined  relationship  between  the  settling  velocity  and  SSC  for  flocculation  and  hindered 
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settling  regions.  Here,  a modified  relationship  of  Hwang,  which  includes  the  effects  of 
flocculation,  salinity  and  temperature  on  settling  velocity,  is  used: 


where  a,  b,  a and  P are  sediment-dependent  empirical  coefficients,  0 is  the  temperature,  v(6) 
is  the  temperature-dependent  fluid  kinematic  viscosity,  p(0,5,c)  is  the  temperature,  salinity 
and  SSC  dependent  fluid  density,  and  F(d)  is  a temperature  fimction  which  reflects  the  effect 
of  temperature  on  flocculation  defined  as  (o^ , where  o)^  is  the  flocculation  settling 

velocity  at  15  °C.  By  reprocessing  the  experimental  data  of  Lau  (1994)  at  different 
temperatures  (see  Appendix  C),  an  empirical  relationship  for  F(d)  is  obtained  as  (Figure  3.2): 


3.4.1  Baroclinic  Effects 

Sediment-induced  stratification  is  considered  in  the  hydrodynamic  equations,  i.e.,  the 
second  and  third  terms  on  the  right  hand  sides  of  Eqs.  (3.2)  and  (3.3),  where  the  bulk  density 
(p)  of  water/ sediment  mixture  are  related  to  SSC  (=c)  by 


(3.30) 


F(0)=1.776-O.O5180,  for  0-O-3O°C 


(3.31) 


3.4  Flow-Sediment  Coupling 


(3.32) 


3.4.2  Vertical  Momentiun  and  Mass  Diffusion  Coefficients 


When  the  flow  is  stratified,  the  buoyancy  effect  tends  to  restore  vertically  moved 
fluid  lumps  back  to  their  original  positions,  and  thereby  causes  a reduction  of  the  turbulent 
transfer  of  momentum  and  mass.  From  the  theoretical  relationship  derived  by  Rossby  and 
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Montgomery  (1935),  Munk  and  Anderson  (1948)  proposed  following  generalized  semi- 
empirical  formulas; 


Figure  3.2.  Dimensionless  median  settling  velocity  as  a function 
of  temperature,  where  o)^  is  the  median  settling  velocity  at  15  °C. 


Momentum  diffusion  coefficient: 


Mass  diffusion  coefficient: 

(3.34) 

where  a^,  b^,  Yp  ^2  Y2  empirical  constants,  and  and  respectively  are 


the  vertical  momentum  and  mass  diffusion  coefficients  in  homogenous  flow.  In  terms  of  the 
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well-known  mixing  length  concept  of  Prandtl  (1925),  and  have  simple  forms  for 

neutral  turbulent  diffusion: 

Momentum  diffusion  coefficient: 


A -I 

^vO  ‘mO 


du 

dz 


=KW^//(-0)(l  +o) 


(3.35) 


Mass  diffusion  coefficient: 


^vO  ^cO^mO 


du 

dz 


=KM,/f(-0)(l  +o) 


(3.36) 


Here  a={z-(,)/H,  1^^  and  are,  respectively,  the  momentum  and  mass  mixing  lengths  in 


a homogeneous,  non-cohesive  flow,  and  Ri  is  the  Richardson  number  defined  as 


Ri^ 


_ g dp/dz 
P (du/dzf- +{dv/dzy 


(3.37) 


Jobson  and  Sayre  (1970)  noted  that  vertical  mixing  of  suspended  sediment  in  open 
channel  flow  occurs  as  a result  of  at  least  two  semi-independent  processes  which  are  shown 
to  be  additive.  These  processes  are:  (1)  Diffusion  due  to  tangential  components  of  turbulent 
velocity  fluctuations,  which  is  the  predominant  turbulent  mixing  process  for  fine  sediment 
particles  in  general,  and  for  all  sediment  particles  in  flows  without  strong  vortex  activity;  and 
(2)  diffusion  due  to  centrifugal  force  arising  from  the  curvature  of  fluid  particle  path  lines, 
which  is  significant  for  coarse  sediment  in  flows  with  strong  vortex  activity.  They  derived 
theoretically  the  following  expressions  for  the  total  turbulent  transfer  coefficient  of  sediment: 
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(3.38) 


in  which  the  first  term  represents  turbulent  transfer  of  sediment  due  to  rectilinear  velocity 
fluctuations,  and  the  second  term  represents  turbulent  transfer  of  sediment  due  to  curvature 
of  fluid  particle  pathlines.  Coefficients  and  are  assumed  to  be  functions  of  the  particle 

characteristics.  Using  the  method  of  average  curve  fitting  technique,  Jobson  and  Sayre  (1970) 
reported  that  p,  was  0.98  and  0.49  and  a was  0.038  and  0.1  for  fine  and  coarse  sediments, 

respeetively. 

A drawback  of  Eqs.  (3.33)  and  (3.34)  is  that  the  turbulent  diffusion  coefficients  will 
become  zero  if  and  where  Richardson  number  tends  to  infinity,  which  is  the  case  when  the 
vertical  gradient  of  velocity  is  zero.  However,  this  is  not  realistic,  since  diffusion  is 
practically  never  equal  to  zero.  Hence,  to  overcome  this  disadvantage,  an  additional  term 
interpreted  a “background”  value  of  the  diffusion  coefficient  can  be  introduced  in  the 
formulas  of  Munk  and  Anderson  (1948),  i.e.. 

Momentum  diffusion  coefficient: 


(3.39) 


Mass  diffusion  coefficient: 


(3.40) 
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where  and  are  the  background  values  of  the  turbulent  diffusion  coefficients  of 
momentum  and  mass,  respectively.  Representative  values  of  a^,  b^,  y,,  a^,  and  are 
given  in  Table  3.1. 


Table  3.1.  Parameters  for  momentum  and  mass  diffusion  coefficients  in  a stratified  flow 


«1 

Y, 

*2 

Y2 

Source 

1 

60-160 

-0.5 

— 

— 

— 

Rossby  and  Montgomery  (1935) 

1 

10 

-0.5 

1 

3.33 

-1.5 

Munk  and  Anderson  (1948) 

1 

10-15 

-1 

— 

— 

— 

Kent  and  Pritchard  (1959) 

3.5  Solution  Techniques 

To  solve  hydrodynamic  Eqs.  (3.1),  (3.2)  and  (3.3)  and  sediment  transport  Eq.  (3.15), 
a semi-implicit  finite  difference  method  is  applied,  which  discretizes  the  convective  and 
diffusive  terms  by  an  Eulerian-Lagrangian  scheme  (Casulli  and  Cheng,  1992).  This  solution 
method  has  the  advantage  of  a minimum  degree  of  implicitness,  good  stability  and 
consistency,  and  high  computational  efficiency  at  a low  computational  cost. 

As  shown  in  Figure  3.3,  a spatial  mesh  consisting  of  rectangular  cells,  totally 
MxN>^Lj-  and  each  of  length  ax  and  Ay  and  height  ao,  is  introduced.  Each  cell  is  numbered 

at  its  center  with  indices  i J and  k.  The  discrete  w- velocity  component  is  then  defined  at  half- 
integer i and  integers  j and  k.  Similarly,  the  v-velocity  component  is  defined  at  integers  i and 
k and  half-integer  j.  The  vertical  velocities,  to  and  w,  are  defined  at  integers  i and  j and  half- 
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integer  k.  Water  surface  elevation  C is  defined  at  integers  / and  j.  The  undisturbed  water 
depth  h{x^)  and  the  total  water  depth  are  specified  at  the  both  m and  v points.  Finally, 
SSC,  denoted  by  c,  and  fluid  density  p are  defined  at  integers  i ,j  and  k. 


i-m  i 1+1/2 

;+l/2 


j-\  12- 

ax 

^:v,  H,h 
•:  C,  w,  c,  p 
(a)  Horizontal  mesh 


A.  ♦ 1/2  A ^ 

1 k 

K * 1/2  ’ 

•:  M,  V,  c,  p 

Cl),  w 


(b)  Vertical  mesh 


Figure  3.3.  Schematic  diagram  of  computational  mesh  and  notation. 


3.5.1  Discreti2ation  of  Hydrodynamic  Equations 

The  general,  semi-implicit,  discretization  of  the  continuity  Eq.  (3.1)  and  momentum 
Eqs.  (3.2)  and  (3.3)  takes  the  following  forms: 
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Differential  continuity  equation: 


/ 

AO 


n + l 

^i*M2j,k 


n*\ 


Ay 


V + l/2. 


V'  1/”^*  ” 

2Z^ 

k=\  k=\ 


+1 

ij-m,k 


(3.41) 


Differential  momentum  equation  in  the  x-direction: 


/j+i 


At 


Cn*\  _rn*\ 

i*\j  ^ij 
AX 


+B 


n 

i^Xnj^k 


• *2,/, J 

f n+\  n+1  \ 

^i*\/2j,k  ^i*mj,k-lj 

1 

I^AO^ 

(3.42) 


Differential  momentum  equation  in  the  y-direction: 


«+i 

^ij^\n,k 


At 


rn*\  _r«  + l 
^i*\j  ^ij 

AX 


+B 


n 

V + i/2,)t 


^\j*\l2,k*\l} 

f « + l /i+l  ] 

^ij*l/2,k*\  ^ij*y2,kj 

1 « + l /j  + 1 \ 

^ij*\/2,k  ^ij^l/2,k-lj 

1 

PAO^ 

(3.43) 


Here  the  subscripts  denote  spatial  positions,  superscripts  denote  time  steps,  At,  B is  the  sum 
of  the  Coriolis,  baroclinic,  horizontal  diffusion  and  Bingham  yield  strength  terms  discretized 


by  the  hydrodynamic  values  at  time  step  n,  and  v"  are  respective  values  of  w and  v at 


time  step  «,  and  subscript  p denotes  water  particle  position  that  is  currently  located  at  u or 
V point.  To  obtain  position  p,  it  is  assumed  that  the  velocity  field  (u,  v)  remains  unchanged 
from  the  previous  time  step  n to  the  present  time  step  n+\. 
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Substituting  Eqs.  (3.42)  and  (3.43)  into  Eq.  (3.41),  a linear,  five-diagonal  system  of 
equations  for  the  water  surface  elevation,  C,  is  obtained.  This  system  is  symmetric  and  strictly 
diagonally  dominant  with  positive  elements  on  the  main  diagonal  and  negative  ones 
elsewhere.  Thus  it  is  positive-definite  and  has  an  unique  solution.  A substantial  part  of  the 
computational  time  is  utilized  in  solving  this  linear  system  (Casulli  and  Cheng,  1992).  In 
practice,  this  system  can  be  solved  efficiently  by  the  pre-conditioned  conjugate  gradient 
method  (see  Appendix  B for  details),  which  is  fast  and  requires  a minimum  amount  of 
computer  memory  (Bertolazzi,  1990).  Then  the  velocity  components,  u and  v are  obtained 
from  Eqs.  (3.42)  and  (3.43). 

3.5.2  Discretization  of  Sediment  Transport  Equation 

Similar  to  the  discretization  of  the  hydrodynamic  equations,  the  semi-implicit 
discretization  of  Eq.  (3.15)  is  given  by 


where  the  notations  are  the  same  as  in  Eqs.  (3.42)  and  (3.43),  D is  the  sum  of  the  discretized 
vertical  settling  and  horizontal  diffusion  terms  using  the  flow  conditions  and  SSC  at  the 
previous  time  step  n. 

3.5.3  Discretization  of  Consolidation  Equation 

The  consolidating  bottom  layer  of  thickness  H'  is  divided  into  sublayers,  with 

each  sublayer  defined  by  a concentration  c and  thickness  H'/L^ . At  each  time  step  At,  the 


sediment  transport  model  (COSED-UF)  provides  the  net  sedimentation  (mass  per  unit  area). 
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qAt,  at  each  grid  point.  Also,  at  each  time  step,  qAt  is  introduced  onto  the  consolidating  layer 
in  following  way  (Figure  3.4);  (1)  When  net  deposition  takes  place,  i.e.,  ^^0,  the  initial 

concentration  of  the  deposited  sediment  is  prescribed  as  Cj-  (typically  Cy.=80  kg  m ; Odd 

and  Cooper,  1989),  and  the  corresponding  thickness  aH'  (=qAt/Cj)  is  added  to  the  top  of  the 

consolidating  layer.  Then  the  consolidating  layer  is  redivided  into  sublayers,  and  new 

initial  concentrations  at  each  computational  element  are  obtained  by  cubic  spline 
interpolation.  (2)  In  the  case  of  net  erosion  an  erosion  depth,  aH'  , is  subtracted  from  the  top 
of  the  consolidating  layer.  Then,  through  redivision  and  cubic  spline  interpolation,  new  initial 
concentrations  at  each  element  are  obtained  as  before.  Following  this  procedure,  the 
consolidation  process  is  modeled  using  the  discretized  form  of  Eq.  (3.22).  To  increase  the 
modeling  accuracy,  a higher  time  resolution  is  applied,  i.e.,  the  consolidation  time  step 
At' = At /N^,  taking 

An  explicit  scheme  is  used  for  descretizating  Eq.  (3.22): 


Here  the  subscripts  denote  spatial  positions,  superscripts  denote  time  steps  and  ao  ' is  the 
vertical  step  length.  Then,  at  each  time  step  the  total  thickness  of  the  consolidating  layer  is 
calculated  from  mass  conservation  according  to: 


(3.45) 


0 


(3.46) 
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Deposition 
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Figure  3.4.  Schematization  of  the  simulated  consolidation  process, 
where  (a)  is  the  original  consolidating  layer,  (b)  is  the  case 
of  net  deposition  and  (c)  is  the  case  of  net  erosion. 


3.5.4  Properties  of  the  Finite-Differential  Equations 

The  accuracy,  stability,  numerical  diftusion  and  spurious  oscillations  of  the  finite- 
differential  equations  (3.41),  (3.42),  (3.43)  and  (3.44)  depend  on  the  discretization  scheme 
of  the  convective  terms,  when  an  Eulerian-Lagrangian  approximation  is  adapted.  Actually, 
as  expressed  in  Eqs  (3.2),  (3.3)  and  (3.15),  the  convective  terms  can  be  rewritten  more 

compactly  as  Lagrangian  derivatives  according  to 

DG  dG  dG  dG  dG 

= +u +V +(0 

Dt  dt  dx  dy  da 


(3.47) 
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where  G denotes  any  physical  quantity,  e.g.,  the  velocity  or  SSC.  Then  the  Eulerian- 
Lagrangian  scheme  discretizes  the  convective  terms  as 


At 


where  subscript  o denotes  the  current  position  of  a water  particle. 

To  obtain  values  of  Gp  , or  v^”  and  Cp  in  Eqs.  (3.42),  (3.43)  and  (3.44),  the 

Eulerian-Lagrangian  scheme  uses  the  back-tracing  approach  incorporated  by  a suitable 
interpolation  method  using  three  or  more  mesh  points  (Casulli  and  Cheng,  1992),  in  which, 
as  stated  before,  it  is  assumed  that  the  velocity  field  («,  v)  remains  unchanged  from  the 
previous  time  step  n to  present  time  step  n+\.  Here,  these  values  are  approximated  by  a 
bilinear  interpolation  over  the  eight  surrounding  mesh  points  (see  Appendix  B for  details). 
The  back-tracing  time  interval  At"=AtlNp,  where  Np  is  the  back-tracing  step  at  each  modeling 
time  step,  is  usually  selected  to  be  ^5.  In  this  event,  the  Eulerian-Lagrangian  scheme 
becomes  free  of  spurious  oscillations.  Moreover,  numerical  diffusion,  which  can  be  regarded 
as  the  interpolation  error,  is  reduced  when  compared  with  numerical  diffusion  induced  by 
the  “up-wind”  method.  Further  reduction  in  artificial  diffusion  can  be  obtained  by  decreasing 
the  spatial  steps  ax.  Ay  and  ao  and  increasing  back-tracing  step,  Np . Complete  elimination 

of  numerical  diffusion  can  be  achieved  by  using  a higher-order  interpolation  formula,  but  the 
resulting  method  may  introduce  some  spurious  oscillations.  Applications  of  this  scheme  to 
problems  with  large  vertical  diffusion,  or  K^,  or  small  vertical  spacings,  ao,  have 


54 


suggested  the  use  of  an  implicit  discretization  only  for  the  vertical  diffusion  terms.  In  fact, 


it  can  be  shown  that  the  stability  condition  for  the  above  scheme  is  simply  given  by 


(Greenspan  and  Casulli,  1988) 


A^,  K^,  Ky 


-1 


(3.49) 


Evidently,  when  A^=Ay=K^-Ky=Q ,ibis  scheme  becomes  unconditionally  stable.  However, 

if  one  restricts  the  back-tracing  process  to  within  one  cell,  i.e.,  the  back-traced  water  particle 
is  not  allowed  to  emerge  out  of  the  cell  where  it  would  start  from  the  beginning,  the  stability 
condition  becomes 


where  and  respectively  are  the  maximum  velocities  in  the  x and>^  directions  within 
the  modeled  domain. 


3.6.1  Hydrodynamics 

In  this  section,  results  of  the  COHYD-UF  model  will  be  compared  with  two 
analytical  solutions.  The  first  comparison  will  primarily  investigate  the  model’s  capability 
to  simulate  the  time  propagation  step.  The  second  comparison  will  attempt  to  validate  the 
model’s  ability  to  compute  nonlinear  effects. 


(3.50) 


3.6  Basic  Simulations 


1 . Comparison  with  a Linear  Analytical  Solution:  Consider  tidal  flow  through  an 
open  channel  connected  to  the  sea  at  the  mouth  (x=0)  and  closed  at  the  uphead  end  (x=/). 
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Neglecting  the  advective  terms,  bottom  friction  and  wind  surface  stress,  the  one-dimensional 
hydrodynamic  equations  for  flow  through  the  channel  are  (Ippen  1966) 

Momentum: 


du  ac  n 

— +g— =0 

dt  dx 


(3.51) 


Continuity: 


dt  dx 


(3.52) 


where  U is  the  vertical  mean  velocity  in  the  x direction  and  depth  h is  assumed  to  be 
constant.  The  selected  boundary  conditions  associated  with  Eqs.  (3.51)  and  (3.52)  are: 

At  the  uphead  end  of  the  basin: 

U{l,t)=Q  (3.53) 

At  the  mouth  of  the  basin: 

C(0,0=v4oSino)or  (3.54) 

where  and  0)q  respectively,  are  the  amplitude  and  angular  frequency  of  the  forcing  tide 

at  the  open  boundary.  Selecting  a uniform  rectangular  cross-section  of  the  channel  and  only 
considering  the  first  mode  of  oscillation,  the  solutions  of  Eqs.  (3.51)  and  (3.52)  are: 

Water  surface  elevation: 


C(V)= 


AqCosHI-x) 

coskl 


sincOpt 


(3.55) 


Velocity: 
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Af.C  sink(/-x) 

cosWq/  (3.56) 

hcoskl 

where  is  the  wave  speed  ( -y/^)  and  k is  the  wave  number  ( =g)q/C^). 

To  test  the  hydrodynamic  model  results  against  the  above  solutions,  a rectangular 
basin  with  a constant  water  depth  of  20  m and  basin  length  of  59  km  is  considered.  Assume 
that  a periodic  tide  with  an  amplitude  of  1 m and  a period  of  7=  12  hr  is  forced  at  the  mouth 
of  the  basin.  The  numerical  solution  is  obtained  by  discretizing  the  basin  into  30  grids  with 
ax=2  km  and  time  step  At=20  min. 

Figure  3.5  shows  the  water  surface  elevation  near  the  mouth  (x=10  km),  at  the  mid- 
point of  the  basin  (x=30  km)  and  near  the  closed  boundary  (x=58  km).  Likewise,  Figure  3.6 
presents  the  velocity  near  the  mouth  (x=l  1 km),  at  the  middle  point  (x=3 1 km)  and  near  the 
closed  boundary  (x=51  km).  Both  Figures  3.5  and  3.6  demonstrate  that  there  is  reasonably 
good  agreement  between  theoretical  and  numerical  solutions,  even  though  a slight  phase  shift 
of  velocity  between  the  theoretical  and  numerical  results  occurs  because  terms  higher  than 
zero  order  are  neglected  in  the  analytical  solution. 

2.  Comparison  with  a Non-linear  Analytical  Solution:  When  nonlinear  terms  are 
included  in  the  one-dimensional  shallow  water  equations,  it  is  not  possible  to  obtain  an  exact 
solution.  However,  one  can  use  harmonic  analysis  to  develop  an  analytical  solution,  which 
is  still  cumbersome  because  the  high  order  terms  are  difficult  to  solve  for.  Accordingly,  here 
only  the  zeroth  and  first  order  harmonic  solutions  are  considered  (Liu,  1988). 
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Without  the  inclusion  of  Coriolis  and  bottom  friction  forces,  the  governing  equations 
for  one-dimensional  nonlinear  tidal  motion  can  be  written  as 


dU  j.dU  dC  „ 
— +U — +^— =^0 

dt  dx  dx 


(3.57) 


Equation  (3.52)  remains  unchanged.  The  boundary  conditions  and  the  geometry  of  the  basin 
considered  are  the  same  as  before.  Neglecting  solution  terms  higher  than  first  order,  the 
solutions  are  (Liu,  1988): 

Water  surface  elevation: 


Af.cosk(l-x) 

C(x,/)=-5 ^ 

coski 


smwr 


A^k 


ihcos^kl 


Velocity: 


xsin2^/  -x) + — - — (sin2A:(/  +x) -tai)2klcos2k{l  -x)) 
cos4kl 


rrr  a 

U{x,t)= ^ COSO)t+- 


(3.58) 


cos2o)/ 


hcoskl 


%h  ^co^kl 


xcoslkQ  -x)+—s\tQ.kil  -x)  - — ^ — [coslkil  +x)  -XzrQ.klsirQ.Jdl  -x)) 
2k  cos4kl 


(3.59) 


sin2o)/ 


The  same  forcing  condition  at  the  basin  mouth  as  before  is  used  to  compare  the 
solutions  with  model  results.  Figure  3.7  shows  the  comparison  of  water  surface  elevation, 
while  Figure  3.8  presents  the  velocity.  It  can  be  see  that  the  model  results  are  reasonably 
close  to  the  analytical  solutions.  As  before,  the  slight  phase  shift  of  velocity  between  the 
theoretical  and  numerical  results  occurs  because  terms  higher  than  first  order  have  been 
neglected  in  the  analytical  solution. 
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Figure  3.7.  Modeling  ID  non-linear  hydrodynamic  equation  Figure  3.8.  Modeling  ID  non-linear  hydrodynamic  equation 

for  tidal  flow  in  an  open  channel.  Lines  are  simulations  and  for  tidal  flow  in  an  open  channel.  Lines  are  simulations  and 

open  circles  represent  analytical  solutions.  open  circles  represent  analytical  solutions. 
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3.6.2  Sediment  Transport 

To  check  the  applicability  of  the  COSED-UF  numerical  scheme,  the  following 
modeling  tests  against  some  special  forms  of  the  governing  equations  having  exact  solutions 
as  well  as  field  data  in  a river  are  carried  out.  Plots  of  numerical  and  analytical  solutions  are 
accordingly  presented. 

1 . Steady  State  One-dimensional  Convection-diffusion:  The  basic  equation  is 


at 

dx 


(3.60) 


where  is  the  constant  source-sink  term  and  / is  the  length  of  the  system. 


The  boundary  conditions  are  prescribed  as 


ax 


=0 


lx=/ 


(3.61) 


The  analytical  solution  is 

u 

where  Np^  is  the  Peclet  number,  ul/K^,  which  is  the  ratio  of  convective  transport  to 
diffusive  transport. 

A rectangular  grid  with  equal  lengths  of  the  spatial  step  was  used  in  the  numerical 
solution.  Values  of  the  parameters  used  were:  q^=5,  l=\,  m=1,  Cq=1  and  A!^q=1.  Figure  3.9 

shows  the  comparison  between  the  modeled  results  and  the  analytical  solution. 


61 


2.  Laplace  Equation:  The  Laplace  equation  was  solved  for  a rectangular  domain, 
0<x<l^,  0<y^ly,  with  a parabolic  boundary  condition  for  a quantity  such  as  temperature, 


specified  on>^0. 

The  basic  equation  is 


K^+K—^0 


(3.63) 


The  boundary  conditions  are  given  by 


K. 


c{x,0)-—x{l  -x),  c-Q  on  other  faces 


By  taking  K^=K  -K^,  the  exact  solution  takes  the  form 


/ ^ 4ii:o(l -cos«7r) 

c{x^,t)=2_^  — smi 


/ N 
mix 


n=i  s\r\h{rml  /l  )n^TZ^ 


y X' 


sinh 


V * / 


rm{ly-y) 


(3.64) 


(3.65) 


The  values  1^-3,  1^=4  and  ^"^=40  were  used  in  the  test  problem.  Figure  3.10  shows 

the  comparison  between  the  modeled  results  and  the  analytical  solution. 

3.  One-dimensional  Transient  Heat  Conduction:  The  basic  equation  is 


5c  5^c 


at  0<x</ 


(3.66) 


The  initial  condition  is 


c=0,  0<x^/,  /=0 


(3.67) 
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Figure  3.9.  Modeling  ID  convection-diffusion  equation. 
Line  is  analytical  solution  and  dots  represent  model  simulations. 
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Figure  3.10.  Modeling  2D  Laplace  equation.  Lines  are  analytical 
solutions  and  open  circles  represent  model  simulations. 
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The  boundary  conditions  are 


c(l,t)=c.,  ^ 


dx 


=0 


x=0 


(3.68) 


The  analytical  solution  is  (Carslaw  and  Jaeger,  1959) 


^ = \-ly  illZe  -(2»-i)VV4^Q.(2n+l):tx 
^ 21 


Cq  7I„=o  2n+l 


(3.69) 


where  is  the  elapsed  time  with  T^i=Kj/l^.  The  values  Cq=1,  1=2  and  A:^=0.5  were  used 

in  the  test  problem.  Comparison  between  the  modeled  results  and  the  analytical  solution  is 
shown  in  Figure  3.11. 

4.  Heat  Conduction  with  Radiation:  This  problem  was  chosen  to  check  the  flux 
boundary  condition  formulation  necessary  in  the  model  due  to  the  resuspension  term  at  the 
bed.  The  governing  equation  and  initial  condition  are  the  same  as  Eqs.  (3.66)  and  (3.67), 
respectively  with  the  boundary  conditions 


c(/,0=0. 


dc 
— +a,c 

dx  ‘ 


=0 


x=0 


(3.70) 


The  exact  solution  is  (Carslaw  and  Jaeger,  1959) 


/ \ 
1 +C£^ 

» 2(p 

sin 

i x] 
1-- 

c _ 

-V 

~ n 

1 1) 

^0 

1 +al 

\ ‘ y 

n = l 

P 

„(a/+a?/2+p^) 

(3.71) 


where  is  the  linear  heat  transfer  coefficient,  is  the  linear  heat  transfer  coefficient, 
is  the  positive  roots  of  Pcotp+a^^O  and  is  the  elapsed  time  as  defined  before.  The 
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Figure  3.1 1.  Modeling  ID  transient  heat  conduction.  Lines 
are  analytical  solutions  and  dots  represent  model  simulations. 
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Figure  3.12.  Modeling  heat  conduction  with  radiation.  Lines 
are  analytical  solutions  and  dots  represent  model  simulations. 
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values  CL -0.7,  c^-l  and  /=1  were  used.  Comparison  between  the  modeled  results  and  the 

analytical  solution  is  shown  in  Figure  3.12. 

5.  Transient  Convection-diffusion:  The  governing  equation  is 


dc  dc  „ d^c  „ A ^ 

— +u — -K^ =0,  at 

dt  dx 


(3.72) 


The  boundary  conditions  are 

c(0,/)=Co  , c(~,/)=0  (3,73) 

By  taking  the  initial  condition  as  c(x,0)=0,  the  analytical  solution  takes  the  form 


(Ogata  and  Banks,  1961) 


c{x,t) 


=exp 


ux  u t 


/erfc 


4K. 


-exp 


u\t-0 


4K 


-t-erfc 


X 


(3.74) 


Here  ^ is  the  integration  variable  and  erfc  is  the  complementary  error  function.  The  values 
u=\  and  .^^=0.5  were  used  in  the  modeling  test.  Comparison  between  the  modeled  result  and 

the  analytical  solution  is  shown  in  Figure  3.13. 

6.  Three-dimensional  Laplace  Equation:  The  three-dimensional  equation  was  solved 
for  a cubic  domain,  0<x</^,  0iy<l^,  0<z</^  with  a power  boundary  condition  for  a quantity 

such  as  temperature,  specified  on  z=0. 

The  governing  equation  is 


d^c 

K —+K  —+K  —=0 
^-.2 


dx" 


dy" 


(3.75) 
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Figure  3.13.  Modeling  ID  transient  conveetion-diffusion  equation. 
Lines  are  analytical  solutions  and  dots  represent  model  simulations. 
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By  taking  K^=K^=K^=Kq,  the  exact  solution  has  the  form 


EE 
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The  values  ly-26,  /^=14  and  /Tq=160  were  used  in  the  modeling  test. 


Comparison  between  the  modeled  results  and  the  analytical  solution  is  shown  in  Figure  3.14. 


Figure  3.14.  Modeling  3D  Laplace  equation.  Lines  are 
analytical  solutions  and  dots  represent  model  simulations. 


7.  Modeling  Test  Against  Field  Observations:  To  further  test  the  sediment  transport 
model,  observed  data  of  SSC  and  velocity  profiles  in  the  Savannah  River  estuary  by 


z//,=2/14  z//,=  3/14  r//„=5/14 
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Ariathurai,  et  al.  (1977)  were  used.  Neglecting  lateral  variations  of  SSC  and  velocity  due  to 
the  fact  that  it  is  a narrow  estuary,  the  model  is  used  in  the  vertical  2D  code.  In  the  modeled 
domain  only  three  grid  points  - located  upstream,  mid-point  and  downstream  of  the  6.7  km 
long  estuarine  segment  considered  - were  set.  The  velocity  profiles  measured  at  these  three 
points  were  inputted  for  using  in  the  sediment  transport  model.  The  SSC  profiles  measured 
at  the  two  end  sites  of  the  segment  were  prescribed  as  boundary  conditions  and  that  at  the 
middle  site  was  used  for  verification.  In  the  simulation,  the  following  effects  were  included: 
longitudinal  and  vertical  turbulent  diffusion,  longitudinal  advection,  sediment  settling, 
suspended  sediment  deposition,  and  bottom  erosion.  The  results  are  shown  in  Figure  3.15. 
It  is  seen  that  the  simulated  results  compare  approximately  with  observations.  It  is  likely  that 
the  degree  of  comparison  may  improve  provided:  (1)  modeling  parameters  are  further  fine- 
tuned,  and  (2)  the  full  3D  code  is  used  along  with  more  accurate  boundary  conditions  for 
SSC. 

3.6.3  Consolidation 

For  testing  the  consolidation  model  two  sets  of  experimental  data  on  consolidation, 
without  and  with  the  net  deposition  flux  at  the  bed-fluid  interface  {q),  are  considered  as 
follows. 

1 • Consolidation  without  Deposition:  Consolidation  without  deposition  at  the  bed- 
fluid  interface  was  examined  by  Toorman  and  Berlamont  (1993)  using  an  estuarine  mud 
from  Doel  Dock  in  Belgium.  In  their  experiments,  which  were  conducted  in  a settling 
column  described  by  Van  den  Bosch  et  al.  (1988;  1989;  1990),  both  the  consolidation  rate 
and  the  density  profiles  of  the  consolidating  layer  were  measured.  The  rate  of  consolidation 
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takes  the  form  of  Eq.  (3.23),  for  which  Toorman  and  Berlamont  (1993)  selected 
m s-‘,  m s-^  c^j=20  kg  m-^  c^2=205  kg  m'^,  c=160 

kg  m'^,  m=3and  n=\3  (Figure  3.16),  along  with  the  initial  concentration  of  deposited 

sediment  Cy=80  kg  m . The  results  are  shown  in  Figure  3.17  for  the  vertical  profiles  of 

sediment  concentration  within  the  consolidating  layer,  and  in  Figure  3.18  for  the  time 
variation  of  surface  elevation  of  the  consolidating  layer,  where  is  the  initial  thickness  of 

the  layer.  Both  Figures  3.17  and  3.18  show  that  the  modeled  results  compare  favorably  with 
the  experimental  data. 

2.  Consolidation  with  Deposition:  Experiments  on  consolidation  with  deposition  at 
the  bed-fluid  interface  were  carried  out  by  Burt  and  Parker  (1984).  In  their  tests,  fixed  masses 
of  an  estuarine  mud  at  a fixed  concentration  were  added  at  24  hour  intervals  to  a settling 
column  of  10  m height  and  0.092  m diameter.  A total  of  7 beds  were  added;  density  profiles 
were  measured  24  hours  after  each  bed  was  added  and  immediately  prior  to  the  next  addition. 
Density  was  also  measured  15  days  after  the  last  bed  was  added.  The  material  added  for  each 
layer  comprised  of  4 liters  of  suspension  at  26.3  kg  m solids  content  giving  a dryweight 


mass  of  105.2  g solids. 


Elevation  above  bottom  (m) 
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Figure  3.15.  Modeling  SSC  (unit:  kg  m Solid  lines  are  simulations  and 
dashed  lines  represent  field  data  observed  from  1800,  9/24/68  to  0400,  9/25/68 
in  Savannah  River  estuary  (after  Ariathurai,  et  al.,  1977).  Contours  from  the 
surface  to  the  bottom  are  0.1,  0.25,  0.5,  1,  1.5  and  2,  respectively. 


-5 


Figure  3.16.  Consolidation  rate,  o)^^,  as  a function  of  dry  density  for  Doel 
Dock  mud  with  c^=80  kg  m'^  (after  Toorman  and  Berlamont,  1993). 
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Figure  3.17.  Modeling  laboratory  data  of  Toorman  and  Berlamont  (1993) 
on  consolidation.  Lines  are  model  simulations  and  points  represent  data. 


Figure  3.18.  Modeled  consolidation  curve  (solid  line)  compared  with 
the  laboratory  data  (open  circles)  of  Toorman  and  Berlamont  (1993). 
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Through  calibration,  the  coefficients  in  Eq.  (3.23)  are  taken  as  = 1 x 1 0'^  m s , 

^ ^ ^ ^^8  m'^,c^2=680  kg  m'^,c,=15  kg  m and n,=15 

(Figure  3.19).  is  taken  as  26.3  kg  m In  the  modeling  test,  the  deposited  materials  were 

introduced  into  the  consolidating  layer  in  the  same  way  as  the  experiments.  The  results  of  the 
vertical  profiles  of  sediment  concentration  within  the  consolidating  layer  are  shown  in  Figure 
3.20.  It  is  seen  that  the  modeled  results  are  in  reasonably  good  agreement  with  the 
experimental  data. 


Dry  density,  c (kg  m'^) 

Figure  3.19.  Consolidation  rate,  as  a function  of  dry  density  for 
the  laboratory  tests  of  Burt  and  Parker  (1984)  with  c/=26.3  kg  m'\ 
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Figure  3.20.  Modeling  laboratory  data  of  Burt  and  Parker  (1984) 
on  consolidation  with  deposition  at  the  bed-fluid  interface. 
Lines  are  model  simulations  and  points  represent  data. 


3.6.4  Interfacial  Entrainment 

Modeled  sediment  entrainment  over  the  lutocline  was  tested  against  experimental 
data  of  Mehta  and  Srinivas  (1993)  using  Eq.  (3.26).  In  these  experiments,  values  of 
^jj=3xl0'^  and  D^=1.6xl0'^  were  reported.  In  their  race-track  flume  experiments,  the 
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initially  introduced  kaolinite  mud  layer,  with  a clear  upper  layer  of  water,  was  entrained  as 
the  upper  layer  velocity  was  increased  gradually.  In  the  modeling  test,  the  measured  velocity 
profiles  were  inputted.  The  same  initial  mud  layer  as  in  the  experiment  was  also  inputted. 
Then  the  mud  layer  was  allowed  to  entrain  with  time  using  Eq.  (3.26).  The  modeled  and 
observed  results  of  the  vertical  profiles  of  SSC  at  different  times  are  shown  in  Figure  3.21. 
The  modeled  results  are  seen  to  be  in  reasonable  agreement  with  the  experimental  data. 


Figure  3.21.  Modeling  laboratory  data  on  entrainment  by  Mehta  and 
Srinivas  (1993).  Lines  are  model  simulations  and  points  represent  data. 


CHAPTER  4 

FIELD  INVESTIGATION  AND  DATA  ANALYSIS 
4.1  Study  Area  Description 

Jiaojiang  estuary  is  located  on  the  east  coast  of  China,  about  200  km  south  of  the 
Yangtze  River  (Figure  4.1).  Jiaojiang  is  200  km  long  and  drains  a basin  of  6,500  kml  The 
estuarine  segment  is  only  about  35  km  in  length.  The  mean  water  depth  below  mean  sea  level 
(MSL)  in  this  segment  is  4-7  m,  and  the  mean  width  is  about  1 .2  km  with  a maximum  of  1 .8 
km.  Seaward  of  its  mouth  the  width  increases  rapidly,  forming  Taizhou  Bay  in  shallow 
coastal  waters.  The  annual  inflow  rate  is  about  6.66x10^  m^,  with  a mean  discharge  of 

210  m^  s . Semi-diurnal  macro-tides  prevail  there.  At  the  mouth  the  mean  tidal  range  is 
about  4 m with  a spring  range  of  6.3  m.  The  depth-mean  peak  tidal  current  can  be  up  to  2.0 
ms"'.  The  tidal  wave  is  strongly  distorted  in  the  estuary  within  a short  distance  from  the 

mouth.  Within  the  estuary  itself,  the  duration  of  ebb  exceeds  flood  by  1-2  hr.  As  a result  the 
tidal  currents  exhibit  a asymmetry;  e.g.,  at  Haimen  the  measured  maximum  flood  and  ebb 
currents  are  2.1  and  1.8  m s , respectively  (Zhou,  1986;  Dong,  et  al.,  1997;  Guan,  et  al., 
1998). 
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Figure  4. 1 . Location  map  of  Jiaojiang  estuary,  China.  Depths  are  in  meters 
below  lowest  astronomical  tide.  Ml  and  M2  are  mooring  sites  (Table  4.1); 

Cl,  C2,  C3  and  C4  are  velocity  measurement  and  SSC  profile  sampling 
stations  (Table  4.1);  C6  is  the  site  of  ASSM  (Table  4.1)  and  T1-T6  are  tide 
stations.  The  region  between  the  double  dotted  lines  is  the  modeled  domain. 


Jiaojiang  is  highly  turbid  with  near-bottom  SSC  often  exceeding  10  kg  m and 
depth-mean  SSC  greater  than  5 kg  m’^  (Li,  et  al.,  1993;  Dong,  et  al.,  1997).  Lutoclines 
along  with  fluid  mud  of  1-3  m thickness  regularly  appear  except  around  flood  slack.  The 
river  sediment  is  mostly  in  the  clay  and  fine  silt  size  range.  The  suspended  sediment  is 
dominated  by  clayey-silt  with  a dispersed  mean  particle  size  of  4-6  pm.  The  fluvial  sediment 
is  estimated  to  be  1.23x10^  tons  yr  , with  a mean  inflow  SSC  of  0.18  kg  m i.e.,  about 
4%  the  SSC  in  the  estuary.  Despite  this  high  load,  sediment  discharge  into  Taizhou  Bay  is 
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only  about  1 % that  of  the  Y angtze  river  (1.2x10*  tons  yr ' ' ) (Yu,  1987).  The  sedimentation 

rate  in  this  estuary  is  up  to  0.2  m yr‘‘,  requiring  regular  maintenance  dredging  for 

navigation  from  the  sea  to  the  port  of  Haimen  (Li,  et  al.,  1992).  Haimen  is  located  near  the 
mouth  of  the  estuary  and  is  an  economic  center  of  this  region.  Thus,  it  is  most  important  to 
maintain  the  depth  of  the  navigation  channel.  For  that  reason,  extensive  engineering  and 
scientific  investigations  have  been  carried  out  in  this  estuary.  These  studies  have  focused  on 
tidal  hydrodynamics,  sediment  transport  and  sedimentation  by  means  of  field  investigations 
and  horizontal  and  vertical  2D  numerical  modeling  (e.g.,  Bi  and  Sun,  1984;  Fu  and  Bi,  1989; 
Li,  et  al.,  1993;  Dong,  et  al.,  1997;  Guan,  et  al.,  1998). 

4.2  Experimental  Plan.  Methods  and  Instruments 

In  order  to  examine  flow  and  sediment  dynamics  in  the  Jiaojiang  estuary,  three  field 
campaigns  were  conducted.  Table  4.1  provides  the  dates,  measurements,  locations,  elevations 
at  which  measurements  were  made  and  apparatuses  used.  Sampling  sites  are  shown  in  Figure 
4.1.  In  terms  of  the  objectives  of  these  investigations,  the  observations  can  be  classified  into 
four  types,  namely:  (1)  fluid  mud  observations  at  sites  Ml  and  M2,  (2)  lutocline  observations 
at  site  C6,  (3)  vertical  profiles  of  SSC,  currents,  temperature  and  salinity  at  sites  Cl , C2,  C3 
and  C4,  and  (4)  tidal  elevations  at  sites  T1-T6. 

4.2.1  Fluid  Mud  Observations 

Fluid  mud  related  observations  were  carried  out  at  site  Ml  during  April  13-23,  1991, 
and  at  M2  during  November  9-25,  1995.  In  both  investigations,  the  monitoring  system 
consisted  of  a 2.7  m high  metal  frame  containing  an  Inter-Ocean  (model  S4)  vector- 
averaging electromagnetic  current  meter  and  six  self-logging  optical  fiber  analite 


Table  4.1.  Summary  of  Jiaojiang  field  campaigns 
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nephelometers.  The  nephelometers  were  fixed  to  the  frame  at  six  elevations  near  the  bottom 
(Table  4.1).  However,  in  the  1991  campaign,  useful  data  were  recorded  at  only  three 
elevations.  These  transducers  recorded  data  at  intervals  of  5 min  in  1991  and  10  min  in  1995. 
The  data  were  bin-averaged  over  interval  of  1 min  with  sampling  at  1 s interval. 

The  Inter-Ocean  current  meter  was  mounted  on  the  frame  at  1.75  m in  1991  and  1.5 
m in  1995,  respectively,  above  the  bottom.  The  meter  logged  velocity  data  sampled  at  0.5  s 
interval  and  bin-averaged  over  1 min.  The  logged  interval  was  5 min  in  1991  and  10  min  in 
1995,  so  that  the  meter  operated  for  only  1 min  every  5 or  10  min.  Practically  calm  water 
surface  prevailed  throughout  the  experiments  with  negligible  waves. 

4.2.2  Lutocline  Observations 

Lutocline  observations  were  conducted  at  site  C6  during  0600-1600  on  November 
15,  1995.  A ship-borne  Acoustic  Suspended  Sediment  Monitor  (ASSM)  (made  by  the 
Shanghai  Acoustics  Laboratory,  Academia  Sinica)  was  used  to  detect  the  lutocline.  The 
ASSM  consisted  of  a 0.5  MHZ  acoustic  transducer/receiver.  The  acoustic  probe  was 
deployed  1-2  m below  the  water  surface.  The  entire  system  was  under  control  of  a PC  for 
synchronization  of  sampling,  preliminary  data  reduction  and  storage.  The  device  had  a pulse 
length  of  about  40  ps,  and  measured  the  vertical  profiles  of  sound  scattered  from  suspended 
sediments  in  the  range  bins  at  0.6  sec  interval  with  a vertical  resolution  of  5 cm.  The  data 
were  sampled  at  a rate  of  approximately  75  kHz  for  9 min  bursts.  Each  data  burst  consisted 
of  900  profiles  of  backscattered  acoustic  energy  from  suspended  sediment  particles  between 
the  bed  and  the  acoustic  probe. 
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4.2.3  Observations  of  SSC.  Currents.  Salinity  and  Temperature 

Vertical  profiles  of  SSC,  tidal  currents,  temperature  and  salinity  were  obtained  at  sites 
Cl,  C2,  C3  and  C4  (Figure  4.1),  respectively  using  Niskin  bottles  (each  60  cm  long  and  10 
cm  in  diameter),  SLC9-1  digital  current  meters  (made  by  the  Institute  of  Marine  Instrument, 
Qingdao),  CTD  probes,  turbidimeters  and  a ‘mud  probe’  which  consisted  of  a CTD 
transducer  equipped  with  an  Analite,  infra-red,  backscattering  nephelometer  (Wolanski,  et 
al.,  1988).  This  profiler  was  able  to  measure  SSC  from  0.03  to  80  kg  m The  time  series 

of  such  profiles  were  collected  over  1 or  2 tidal  cycles.  Water  samples  were  collected  in  the 
Niskin  bottles  at  six  elevations  between  the  water  surface  to  the  bottom.  Additionally,  water 
samples  were  collected  at  0.3,  0.6  and  1 m elevations  above  the  bed  using  three  horizontally 
deployed  Niskin  bottles  mounted  on  a solid  frame.  These  bottles  were  raised  on  board  within 
about  half  a minute  of  sampling  underwater,  and  water  samples  were  then  drawn 
immediately  for  analysis  of  SSC,  salinity  and  sediment  size  (Li,  et  al.,  1993;  Dong  et  al., 
1997). 

4.2.4  Tidal  Elevations 

Tides  were  obtained  at  six  sites  from  T1  to  T6  (Figure  4.1)  for  a one-month  period. 
Each  time-series  was  processed  by  harmonic  analysis  (Dong,  et  al.,  1997;  Guan,  et  al.,  1998). 

4.3  Experimental  Data 

4.3.1  Sediment  Size 

The  river  sediment  is  mostly  in  the  clay  and  fine  silt  size  range.  The  suspended 
sediment  is  dominated  by  clayey-silt  with  a dispersed  mean  particle  size  of  4-6  pm  (Figure 
4.2).  About  50%  by  volume  of  the  particles  constitute  very  fine  silt  (size  range  between  4-16 
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urn).  The  particle  fraction  with  size  less  than  4 pm  accounts  for  about  40%  of  the  material 
by  volume  (Li,  et  al.,  1993;  Li,  et  al.,  1999). 
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Figure  4.2.  A representative  frequency  distribution  of 
suspended  sediment  (dispersed)  size  in  the  Jiaojiang  estuary. 


4.3.2  Tides 

1.  Elevation:  Harmonic  analysis  showed  that  the  Mj  and  S2  constituents  dominate 

this  estuary.  Figure  4.3  shows  the  time  series  of  tidal  elevations  at  sites  T1  and  T5  during  a 
spring  tide.  It  is  seen  that  the  tide  was  distorted  as  it  propagated  upstream  from  the  mouth. 
For  example,  at  T1  (outside  the  mouth),  the  flood  lasted  5.78  hr  and  the  ebb  6.64  hr.  In 
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contrast,  at  T5  (inside  the  estuary)  the  flood  was  5.29  hr  and  ebb  was  7. 13  hr.  In  other  words, 
within  a short  distance  of  about  1 8 km  from  the  mouth  to  upstream,  the  duration  of  flood 
decreased  by  about  one-half  hour. 


Figure  4.3.  Time  series  of  tidal  elevation  at  sites  T1  and  T5  during 
a spring  tide  from  0000  hr  on  1 1/05/94  to  0100  hr  on  1 1/06/94. 


2.  Currents:  As  a result  of  tidal  distortion,  the  tidal  currents  also  exhibited  an 
asymmetry  as  shown  in  Figures  4.4-4.6.  It  is  seen  that  the  flood  current  duration  was  shorter, 
and  its  strength  was  greater,  than  the  corresponding  ebb  current  values.  Thus,  for  example, 
at  C4  the  ratio  of  depth-mean  maximum  flood  current  to  ebb  current  was  about  1 .4  during 
the  spring  tide  and  about  1 .25  during  the  neap  tide. 
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4.3.3  Profiles  ofSSC 

Figures  4.7  and  4.8  show  typical  time  series  of  SSC  profiles  at  C4  during  a spring  and 
a neap  tide,  respectively,  and  Figure  4.9  for  C6  during  a neap  tide.  It  is  seen  that  during  low 
water  slack  and  neap  tide,  the  lutocline  was  distinct  and  well  defined.  In  contrast,  during 
spring  tide  and  peak  flows,  the  lutocline  became  indistinct  even  though  the  fluid  mud  layer 
(as  defined  in  Section  2.1)  became  thicker.  This  observation  supports  the  conclusion  stated 
in  later  Section  5.2  regarding  the  relationship  between  the  robustness  of  the  lutocline  and  the 
flow  conditions.  Also  at  these  sites,  the  vertical  mean  SSC  was  always  greater  than  5 
kg  m and  often  exceeded  10  kg  m . 

4.3.4  ASSM  Data 

Figures  4.10  (a)  and  (b)  show  two  typical  ASSM  outputs,  each  over  a 60  s time 
segment.  Observe  that  most  of  the  time  the  ASSM  signal  had  a distinct  step  structure.  As 
demonstrated  in  Figure  4.1 1,  this  step  structure  matches  the  real  lutocline  position  detected 
by  the  turbidimeter.  Hence  it  is  reasonable  to  treat  the  step  structure  in  the  ASSM  signal  as 
the  location  of  lutocline.  Also  observed  are  lutocline  undulations,  suggesting  noteworthy 
internal  wave  generation  and  likely  sediment  entrainment  activity.  A smoother  interface  was 
observed  during  ebb  [Figure  4.10(b)]  than  during  flood  [Figure  4.10(a)]. 

The  SSC(=c)  can  be  calculated  from  the  ASSM  reading,  F^,  from  (Thome  et  al., 

1994) 

(4.1) 

where  h'  is  the  depth  below  the  acoustic  probe  and  k^,  h^,  and  are  sediment- 
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Figure  4.4.  Time  series  of  velocity  at  site  C4  during  a spring  tide.  Observations  began  at 
1700  hr  on  November  5,  1994.  Positive  numbers  signify  flood  and  negative  are  for  ebb. 


Figure  4.5.  Time  series  of  velocity  at  site  C4  during  a neap  tide.  Observations  began  at 
0900  hr  on  Novermber  10,  1994.  Positive  numbers  signify  flood  and  negative  are  for  ebb. 
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Figure  4.6.  Time  series  of  velocity  at  site  C6  during  a neap  tide.  Observations  began  at 
0630  hr  on  November  15,  1995.  Positive  numbers  signifie  flood  and  negative  are  for  ebb. 


Figure  4.7.  Time  series  of  SSC  at  site  C4  during  a spring  tide.  Observations  began 
at  1700  hr  on  November  5,  1994.  Shaded  area  includes  SSC  greater  than  20  kg  m'^. 
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Figure  4.8.  Time  series  of  SSC  at  site  C4  during  a neap  tide.  Observations  began 
at  0900  hr  on  November  10, 1994.  Shaded  area  includes  SSC  greater  than  20  kg  m ^ 


8 10  12  14  16  18 

Time  (hr) 

Figure  4.9.  Time  series  of  SSC  at  site  C6  during  a neap  tide.  Observations  began 
at  0600  hr  on  November  15,  1995.  Shaded  area  includes  SSC  greater  than  20  kg  m 'I 
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(b) 

Figure  4.10.  Typical  raw  ASSM  records  during  the  neap  tide  on  Nov.  15, 
1995,  with  a horizontal  time  scale  of  1 min  and  a vertical  distance  scale 
of  1.25  m.  (a)  was  observed  during  a flood  with  a value  of  Richardson 
number  Ri^  of  about  2,  and  (b)  during  an  ebb  with  Ri^  of  about  150. 
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dependent  constants.  As  mentioned  above,  the  location  of  the  largest  vertical  gradient  of 
can  be  taken  as  the  lutocline  elevation,  . above  the  bottom.  From  the  time  series  of  C^(0 , 

a spectral  analysis  of  internal  waves  can  be  carried  out,  and  related  to  local  hydrodynamics. 


Figure  4.1 1.  Relationship  between  lutocline  elevations  above  bottom  detected  by  the 
turbidimeter  and  by  the  ASSM.  Data  were  collected  during  0600-1600,  Nov.  15,  1995. 


Figure  4.12  shows  three  typical  time-series,  i.e.,  examples  (a),  (b)  and  (c),  of 
lutocline  elevation,  C,Jj),  traced  from  9 min  long  ASSM  records.  It  is  evident  that  there  are 

two  types  of  internal  waves  riding  on  the  lutocline,  namely,  low  frequency  internal  waves 
with  a period  on  the  order  of  1 min  and  high  frequency  internal  waves  with  a period  of  about 
5 s. 
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Figure  4.12.  Time  series  of  lutocline  elevation  at  site  C6  during  a neap  tide  using  ASSM  on  November 
15,  1995.  (a)  and  (b)  were  sampled  during  flood  with  a value  of  Ri^  of  about  1,  and  (c)  during  ebb 
with  Ri^of  about  150.  Solid  lines  are  instantanious  elevations  and  dashed  lines  are  mean  trends. 
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Two  distinct  types  of  internal  waves  on  the  lutocline  have  also  been  reported  by  other 
researchers.  From  the  observation  of  internal  waves  on  the  lutocline  in  a race-track  flume 
experiment,  Mehta  and  Srinivas  (1993)  described  two  types  of  internal  waves;  One  was 
characterized  by  high  frequency  and  breaking.  While  the  other  included  large-amplitude 
solitary  waves  which  interspersed  with  the  breaking  waves  and  seemed  to  decay  without 
breaking.  Shi  (1998)  also  observed  two  types  of  internal  waves  on  the  interface  between  the 
mobile  and  stationary  fluid  mud  in  Hangzhou  Bay  in  China  using  an  ASSM  meter.  These 
included  low  frequency  waves  with  a period  on  the  order  of  minutes  and  high  frequency 
waves  with  a period  of  seconds.  Wright,  et  al.  (1988)  observed  low  frequency  (-6x10"'* 
rad  s‘*)  and  relatively  large  amplitude  internal  waves  on  the  upper  interface  of  the 
underflows  over  the  active  delta  front  of  the  Yellow  River  in  China.  They  found  that  these 
internal  waves  had  frequencies  near  the  local  Brunt-Vaisala  frequency  (~  5 .3  x 1 0 rad  s ‘ ‘ ). 

For  the  purpose  of  spectral  analysis  of  the  high  frequency  internal  waves,  the 
lutocline  elevation  time-series  processed  for  trend  removal  using  an  approach  used  by  Costa 
(1989),  i.e.,  filtering  out  the  low  frequency  internal  waves.  The  results  after  trend  removal 
are  shown  in  Figure  4.13.  Figure  4.14  shows  an  enlarged  portion  of  the  time-series 
accentuating  the  internal  waves.  It  is  seen  that  the  internal  waves  usually  have  sharp  crests 
and  flat  troughs,  which  is  most  similar  to  the  Holmboe-type  interfacial  instability  (Broward 
and  Wang  1972).  In  experiments  of  Srinivas  (1989),  in  general  both  Kelvin-Helmoltz  and 
Holmboe  modes  of  interfacial  instabilities  appeared.  However,  here  the  Kelvin-Helmoltz 
instability  is  difficult  to  identify.  Due  to  the  interfacial  instability,  a net  upward-asymmetric 
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during  a neap  tide,  where  (a),  (b)  and  (c)  correspond  to  Figure  4.12. 
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mass  mixing  can  be  expected  in  this  case.  Mixing  can  be  further  enhanced  if  the  internal 
waves  begin  to  break. 


Figure  4.14.  Typical  profiles  of  internal  waves  exhibiting  sharp 
crests  and  flat  troughs.  Data  were  taken  from  example  (a)  in 
Figure  4.13.  Wave  heights  range  from  0.07  m to  0.23  m. 


Based  on  the  Wiener-Khintchine  theorem,  the  internal  waves  spectrum,  S'Cto '),  can 
be  calculated  using  (Bendat  and  Piersol,  1971) 


5((0^)=— J R(x)cos((O^T)dv 
^ 0 


with  the  auto-correlation  fimction,  /?(t)  , defined  as 


(4.2) 


1 ^ 

R(x)=lim^  f HJ,t)HJJ^x)dt 
/ J 


(4.3) 


-T 


where  is  the  internal  wave  profile,  x is  the  time  interval  or  shifting  time,  T'is  the 


integration  time-limit  and  o)^  is  the  angular  frequency  (Ochi,  1990).  The  auto-correlation 
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functions  for  examples  (a),  (b)  and  (c)  are  shown  in  Figure  4.15.  The  resulting  respective 
spectral  densities  are  plotted  in  Figures  4.16-4.18.  It  is  seen  that  the  maximum  spectral 
density  is  located  around  the  modal  frequency  o)^=T3  rad  s ''  during  flood,  with  Ri^  of 

about  1 (Figures  4.16  and  4.17),  and  around  w^=l.l  rad  s ''  during  ebb,  with  Ri^  of  about 
150  (Figure  4.18). 

A further  examination  of  internal  waves  was  made  by  relating  the  root  mean  square 
(rms)  of  the  height  of  the  internal  wave  as  well  as  its  modal  frequency  (both  for  high  and  low 
frequency  waves)  to  the  global  Richardson  number  and  the  tidal  range.  The  rms  of  the  high 
frequency  wave  height  and  the  corresponding  angular  frequency  were  taken  over  1 min 
segments  of  the  ASSM  output.  One  segment  was  selected  from  each  9 min  ASSM  burst 
(totally  21  bursts)  for  this  purpose.  The  rms  of  the  low  frequency  wave  height  and  the 
corresponding  angular  frequency  were  taken  over  each  9 min  segment  of  ASSM  burst.  The 
global  Richardson  number,  Ri^,  was  calculated  from  Eq.  (3.26)  using  the  current  velocity 

and  SSC  profiles.  The  lutocline  was  considered  to  be  at  the  near  bottom  elevation  having  the 
maximum  vertical  gradient  of  SSC.  The  average  value  of  Ri^  over  1 min  was  taken  for  the 

high  frequency  waves  and  9 min  for  the  low  frequency  waves.  The  rms  wave  height,  is 

calculated  from 


\ 


1 


2 

an 


(4.4) 


where  H is  the  n wave  height  and  N is  the  total  number  of  waves  over  the  selected  data 


Spectral  density,  S (m  s)  Autocorelation  function,  /?(x)  (m 
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Figure  4.15.  Auto-correlation  function  against  time  interval, 
where  (a),  (b)  and  (c)  correspond  to  Figure  4.12. 
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Figure  4.16.  Internal  wave  spectrum  corresponding  to  example  (a)  in  Figure  4.13. 
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Figure  4.18.  Internal  wave  spectrum  corresponding  to  example  (c)  in  Figure  4.13. 
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segment,  i.e.,  1 min  or  9 min.  The  modal  frequeney  is  caleulated  as 

^a=^r-  (4.5) 

where  is  the  duration  of  each  segment. 

The  processed  results  are  shown  in  Figures  4.19-4.26  and  summarized  in  Table  4.2. 
It  is  seen  that  //  and  both  for  high  and  low  frequencies,  have  reasonably  good 

relationships  Avith  the  global  Richardson  number.  In  general,  the  lower  the  global  Richardson 
number,  the  higher  the  and  values.  This  behavior  is  examined  in  Section  4.4. 

There  appears  to  be  no  identifiable  relationships  between  // , (o^  and  the  tidal  range. 

However,  the  Richardson  number  was  typically  much  higher  during  ebb  than  during  flood, 
because  the  ebb  exhibited  comparatively  more  uniform  velocity  profiles  (Figure  4.6).  As 
observed  in  Figures  4.21, 4.22, 4.25  and  4.26,  internal  waves  had  lower  heights  and  angular 
frequencies  during  ebb  than  during  flood.  During  ebb  the  majority  of  the  data  points  are 
located  below  the  mean  trend-line.  In  contrast,  the  majority  of  the  points  are  located  above 
this  line  during  flood. 

Following  Wright,  et  al.  (1988),  a possible  relationship  between  the  internal  wave 
frequency  and  the  Brunt-Vaisala  frequency  is  examined  here.  The  Brunt-Vaisala  frequency, 
0)^,,  is  defined  as 


p dz 


(4.6) 


Modal  frequency,  co„  (rad  s ‘)  nns  internal  wave  height,  (m) 
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Global  Richardson  number,  Ri^ 

Figure  4.19.  rms  of  high-frequency  internal  wave  height 
as  a function  of  global  Richardson  number. 


Global  Richardson  number,  Ri^ 

Figure  4.20.  Modal  frequency  of  high-frequency  internal 
waves  as  a function  of  global  Rechardson  number. 
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Figure  4.21.  rms  of  high-frequency  internal  wave  height  during  ebb  and  flood. 
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Figure  4.22.  Modal  frequency  of  high-frequency  internal  waves  during  ebb  and  flood. 
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Figure  4.23.  rms  of  low-frequency  internal  wave  height 
as  a function  of  global  Richardson  number. 
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Figure  4.24.  Modal  frequency  of  low-frequency  internal 
waves  as  a function  of  global  Richardson  number. 
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Figure  4.25.  rsm  of  low-frequency  internal  wave  height  during  ebb  and  flood. 
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Figure  4.26.  Modal  frequency  of  low-frequency  internal  waves  during  ebb  and  flood. 
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Table  4.2.  rsm  height  and  angular  frequency  of  internal  waves 
as  functions  of  global  Richardson  number  and  tidal  range. 


Type  of 
internal 
wave 

Parameter 

Relationships 

H^{m) 

M s'*) 

/?/q 

i/^=0.25  -0.041og(R/Q) 

a)^=1.4-0.11og(R/o) 

High 

frequency 

Tidal 

range 

Flood:  9 out  of  10  points  are 
above  the  mean  value  of 
// =0.205  m 

Ebb:  5 out  of  6 points  are 
below  the  mean  value  of 
//j  =0.205  m 

Flood:  6 out  of  9 points  are 
above  the  mean  value  of 
o)^=1.326  rad  s '* 

Ebb:  5 out  of  6 points  are 
below  the  mean  value  of 
(0^=1.326  rad  s '* 

R/q 

7/^=0.4375  -0.05251og(/?/o) 

co^=0.1025-0.01251og(R/o) 

Low 

frequency 

Tidal 

range 

Flood:  1 1 out  of  12  points  are 
above  the  mean  value  of 
//^  =0.382  m 

Ebb:  7 out  of  8 points  are 
below  the  mean  value  of 
// =0.382  m 

Flood:  1 lout  of  12  points 
are  above  the  mean  value 
of  0)^ =0.089  rad  s 

Ebb:  6 out  of  8 points  are 
below  the  mean  value  of 
(0^=0.089  rad  s “* 

Similar  to  the  global  Richardson  number  formula  (3.26),  one  can  approximate  O)  as 


0)  = 
V 


SCJC,-C,) 


N \ Kt 


g' 


(4.7) 


where  g'  is  the  reduced  gravity.  To  estimate  using  Eq.  (4.7),  the  SSC  profiles  (Figure  4.9) 
observed  at  the  same  time  as  the  ASSM  data  were  employed.  The  results  indicate  that  o)  has 


a range  of  0.12-0.26  rad  s' '.  As  shown  in  Figures  4.24,  4.26  and  Table  4.2,  the  low 


frequency  waves  had  frequencies  ranging  from  0.06  to  0.12  rad  s'*  and  a mean  value  of 
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about  0.09  rad  s * . Thus,  the  low  frequency  waves  had  modal  frequencies  near  the  local 
Brunt- Vaisala  frequency.  This  type  of  internal  wave  at  haloclines  or  thermoclines  has  been 
extensively  studied  by  previous  researchers  (e.g.,  Lamb,  1945;  Neumann  and  Pierson,  1966; 
Phillips,  1977).  Theoretically,  in  a stratified  fluid  a vertically-moved  fluid  particle  due  to  any 
disturbance  will  be  restored  to  its  original  level  under  buoyancy  force  and  overshift  inertially. 
Consequently,  it  will  oscillate  about  the  original  level  at  the  Brunt-Vaisala  frequency.  Thus, 

it  is  believed  that  the  low  frequency  wave  detected  in  the  Jiaojiang  exhibited  a natural 
oscillation  about  the  lutocline. 

4.3.5  Floe  Settling  Velocity 

In  order  to  examine  the  characteristic  dependence  of  the  settling  velocity  as  a function 
of  SSC  in  this  estuary,  the  SSC  profile  data  shown  in  Figures  4. 7-4. 9 were  used.  One  may 
conveniently  select  the  SSC  profiles  during  slack  water  for  this  purpose,  since  the  quiescent 
water  assumption  is  approximately  satisfied  then.  From  the  sediment  conservation  Eq.  (3.1 5), 
the  simplified  vertical  equation  for  SSC  transport  is 


dc  1 d(o,c 
dt  H da 


(4.8) 


For  solving  Eq.  (4.8)  to  obtain  the  SSC  profile  c(o,f),  the  initial  and  boundary  conditions  are 
specified  as  follows: 

Initial  condition: 


C(O,0)=Cg(O) 


(4.9) 


Zero  settling  flux  boundary  condition  at  the  water  surface  (o=0): 
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(4.10) 


Zero  settling  flux  boundary  condition  at  the  bottom  (o=-l): 

(4.11) 

Based  on  the  above  boundary  value  formulation,  one  can  approximately  calculate  the 
settling  velocity,  o)^,  as  a function  of  SSC.  This  analysis  has  been  programmed  by  Mehta  and 


Li  (1997).  Accordingly,  the  constants  a,  6,  a and  P in  the  settling  velocity  formula  (3.30)  can 
be  determined  by  best-fitting  Eq.  (3.30)  against  the  data  points  (for  settling  velocity  against 
SSC).  To  determine  the  coefficients  in  settling  velocity  formula  (3.30),  the  effects  of 
temperature  and  salinity  on  flocculation  are  neglected  and  the  settling  velocity  formula  is 
simplified  as 


CO. 


ac 


{c^+b^f 


(4.12) 


Through  best-fitting  of  the  calculated  data  points,  the  values  of  the  coefficients  in  Eq. 

(4.12)  are  obtained  as:  a=0.045,  6=6,  a=1.5  and  P=1.51  during  neap  tide  (Figure  4.27),  and 
0=0.23,  6=10,  a=1.5  and  P=1.8  during  spring  tide  (Figure  4.28).  These  values  are  consistent 
with  c (=SSC)  measured  in  kg  m and  co^  measured  in  m s . It  is  observed  that  the 


maximum  SSC  for  flocculation  settling,  (approximately  equal  to  6),  increases  during 

spring  tide  compared  with  neap  tide,  with  C2=6.0  kg  m during  neap  and  C2=8.4  kg  m 

during  spring.  In  contrast,  the  settling  velocity  during  spring  tide  is  less  than  during  neap  tide 
(Figures  4.27  and  4.28).  Also  shown  in  Figures  4.27  and  4.28  is  the  data  point  of  Li,  et  al. 
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Figure  4.27.  Settling  velocity  as  a function  of  SSC  during  a neap 
tide  from  0900  hr  on  1 1/10/94  to  1000  hr  on  1 1/1 1/94.  Solid  line 
is  the  best-fit  of  the  calculated  data  points  using  Eq.  (4.12). 


Figure  4.28.  Settling  velocity  as  a function  of  SSC  during  a spring 
tide  from  1700  hr  on  1 1/05  /94  to  1800  hr  on  1 1/06/94.  Solid  line 
is  best-fit  of  the  calculated  data  points  using  Eq.  (4.12). 
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(1993)  determined  using  the  Postma’s  ‘pipette’  method  (McCave,  1979)  during  the  1991 
campaign  in  the  Jiaojiang  (Table  4.1).  Finally  note  that  in  both  plots  is  the  free  settling 

velocity. 

As  described  by  Burt  (1984),  the  tidal  range,  or  associated  turbulence  has  two  effects 
on  flocculation  and,  consequently,  the  settling  velocity.  Increasing  turbulence  may  enhance 
flocculation  and  at  the  same  time  limit  the  size  of  floes  that  can  be  sustained.  In  other  words, 
depending  on  its  cumulative  effects  in  enhancing  flocculation  and  limiting  floe  size, 
increasing  turbulence  may  either  increase  or  decrease  the  settling  velocity.  In  the  Jiaojiang, 
limiting  floe  size  apparently  dominated  during  spring  tide,  while  enhanced  flocculation 
occurred  during  neap. 

4.3.6  Erosion  Rate  Constant 

The  erosion  rate  constant,  M^,  in  Eq.  (3.29)  can  be  determined  from  the  vertical 

profiles  of  current  velocity  and  SSC  shown  in  Figures  4.4-4.9.  The  vertical  profiles  chosen 
for  this  purpose  are  those  corresponding  to  periods  when  the  suspended  sediment  mass  in  the 
entire  water  column  increased  gradually,  usually  3-4  hr  during  each  period  of  measurement. 
Also,  the  following  assumptions  are  made: 

(1)  SSC  is  approximately  uniform  along  the  estuary,  so  that  the  longitudinal 
advection  term  in  the  sediment  conservation  equation  (3.15)  can  be  omitted,  i.e.,  udc/dx~0. 

(2)  Within  a short  erosion  period  of  3-4  hr,  the  bottom  shear  strength  can  be 


approximately  taken  as  constant. 
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(3)  The  effect  of  temperature  on  erosion  rate  is  negligible. 

From  the  data  two  parameters  are  obtained,  i.e.,  the  terms  in  Eq.  (3.29)  including  the 

excess  bottom  shear  stress,  exp|-XT^)  and  the  erosion  rate,  m^.  The  bottom  shear 


stress  is  calculated  from  the  vertical  mean  velocity  according  to 

_ 9grif  t/ 

Z/l/3 


(4.13) 


where  rij.  is  the  Manning’s  bed  resistance  coefficient  taken  as  0.015  (Dong  et  al.,  1997),  and 

U is  the  vertical  mean  velocity.  The  bottom  shear  strength  is  taken  as  the  bottom  shear  stress 
at  the  beginning  of  each  selected  time  period.  Also,  according  to  Lee  and  Mehta  (1994),  the 
values  of  x=8  and  A=0.5  are  applicable.  Finally,  the  erosion  rate  is  calculated  from  the  SSC 
profiles  according  to 


m 


e 


At 


(4.14) 


where  C is  the  vertical  mean  SSC,  superscript  n denotes  time  and  At  is  the  time  interval 
between  two  profiles. 

Through  best-fitting  of  the  data  points  (Figure  4.29),  the  erosion  rate  constant, 

0.29  kg  N s ■'  is  obtained.  Vinzon  (1998)  obtained  erosion  rate  constants  by  analyzing 
prototype  data  collected  by  Kineke  (1993)  in  the  same  way.  Her  data  yielded  the  range  of 
M to  be  0.25-0.34  kg  N-‘  s-‘. 

max  ® 
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4.4  Properties  of  Internal  Waves 

The  results  in  Section  4.3.4  show  that  both  the  rms  height,  // , and  the  modal 

frequency,  (o^,  of  internal  waves  decrease  with  increasing  Ri^.  An  attempt  is  made  here  to 

explain  this  behavior  by  referring  to  works  of  previous  researchers.  Also  examined  in  this 
section  are  the  celerity  and  length  of  internal  waves. 

4.4. 1 Effect  of  RIq  on 

Since  high  Ri^  implies  high  buoyancy-induced  stabilization  of  the  lutocline, 
increasing  Ri^  should  correlate  with  decreasing  wave  height.  As  described  in  Section  1 .2, 


this  phenomenon  has  been  observed  by  previous  investigators  in  laboratory  experiments  on 
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lutoclines  (e.g.,  Mehta  and  Srinivas,  1993)  and  on  other  pycnocline  (e.g.,  Chou,  1975; 
Narimousa  and  Fernando  1987).  By  assuming  that  the  interfacial  undulations  are  due  to  the 
energy-containing,  mixed-layer  eddies  impinging  on  the  density  interface,  Narimousa  and 
Fernando  (1987)  established  an  empirical  relationship  between  the  wave  height,  // , and 

as  follows 

-1/2 

Y~  (4.15) 

'%ix 

In  Eq.  (4.15),  is  the  mixed-layer  thickness,  i.e.,  the  upper  layer  depth.  According  to  this 

expression,  decreases  with  increasing  Ri^  with  a slope  of  0.5  on  a log-log  plot.  As  noted 

in  Section  4.3.4,  the  observations  in  Jiaojiang  do  show  that  decreases  with  increasing 

RIq,  although  slopes  much  smaller  than  0.5  were  found  (Figures  4.19  and  4.23,  and  Table 

4.2).  This  difference  in  slopes  between  the  Jiaojiang  and  the  laboratory  results  of  Narimousa 
and  Fernando  (1987)  is  believed  to  be  mainly  due  to  different  physical  scales  and  associated 
hydrodynamic  effects  including  the  degree  of  turbulence  and  eddy  lengths. 

4.4.2  Effect  of  Rif.  on  (o 

In  order  to  examine  the  influence  of  the  Richardson  number  on  the  modal  frequency 
of  internal  waves,  the  work  of  Lamb  (1945)  is  introduced  here.  Lamb  analytically  examined 
internal  waves  at  the  interface  of  two  inviscid  fluids  of  densities  p,  and  beneath  the 

other,  and  moving  parallel  to  the  x-axis  with  velocities  f/,  and  , respectively  (Figure 
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4.30).  By  assuming  both  fluids  to  be  of  unlimited  depth  and  taking  the  wave  profile  as 


i(u)J-hc) 


(4.16) 


Lamb  derived  the  following  expression  of  the  dispersion  relationship  for  the  waves: 


P|^+P2^ 

P1+P2 


± 


g(p2~Pl) 

^Pl+P2) 


P1P2 

(Pl+P2)^ 


1 

2 


(4.17) 


where  k is  the  wave  number.  The  first  term  on  the  right-hand  side  of  Eq.  (4.17)  is  referred 
to  as  the  vertically-averaged  velocity,  U,  of  the  two  layers.  It  is  seen  that  the  values  of 

given  by  (4. 1 7)  are  imaginary  if 


g(p2~Pi)  ^ P2  „1 

P,+P2~2 


and  it  is  also  recognized  that  for  two  fluids  of  nearly  equal  densities,  such  as  water  and  fluid 
mud,  p2/(p,+p2)==0.5. 

It  is  evident  that  under  the  condition  imposed  by  (4. 1 8),  two  possible  cases  can  arise 
with  respect  to  the  sign  of  the  second  term  of  Eq.  (4.17).  Considering  Eqs.  (4.16)  and  (4.17), 
it  is  seen  that  taking  the  plus  sign  the  wave  height  will  dissipate  with  time.  This  inherently 
implies  that  the  interface  will  be  stable.  On  the  other  hand,  if  the  minus  sign  is  taken,  the 
wave  height  will  grow  with  time.  In  other  words,  the  interface  will  be  unstable.  For  the 
present  analysis,  only  the  unstable  mode  is  of  interest,  i.e.,  with  the  minus  sign  relative  to  the 
second  term  in  Eq.  (4.15).  Thus,  as  soon  as  (4.18)  is  satisfied  the  interface  will  become 
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unstable.  This  is  known  as  the  Kelvin-Helmholtz  instability  (Delisi  and  Corcos,  1973).  If 
now  one  considers  the  internal  waves  of  all  likely  wavelengths  in  the  estuarine  setting  such 
as  in  the  Jiaojiang,  it  can  be  concluded  that  sufficiently  short  waves  will  be  present  to  cause 
interfacial  instability.  Therefore,  a two-layered  estuarine  shear  flow  characteristically 
unstable. 


Figure  4.30.  Definition  sketch  of  two-layered  flow  system. 


Based  on  the  above,  holding  p,  and  P2  constant,  Eq.  (4.17)  can  be  restated  as 


(0 


‘u,-vr 


a a eft 


V ) 


(4.19) 


where  /l^-g(p2-p,)/^(p,+P2),  5^=p,P2/(Pi+p2)^  and  //^  is  the  effective  water  depth 


defined  as  the  thickness  affected  by  internal  waves.  Further  holding  the  wave  number  k and 
the  mean  velocity  U constant,  one  may  quantitatively  evaluate  the  influence  of  the  velocity 


Ill 


gradient,  \U^-U^\/H^j^,  on  o)^.  Thus,  Eq.  (4.19)  can  be  expressed  in  terms  of  a stream 
Richardson  number,  Ri  , as 


. a eff 


Ri. 


i_ 

2 


where  5^'=5^(p2-p,)/ Pj  and  Ri^  can  be  conveniently  defined  as 


(4.20) 


P,-  _^P2-Pi)^.,r 


(4.21) 


Thus  Ri^  is  conceptually  analogous  to  Ri^  [Eq.  (3.26)].  From  Eq.  (4.20),  it  is  seen  that 
decreases  as  Ri^  increases,  an  observation  that  is  consistent  with  the  data  in  Figures  4.20  and 


4.24. 

4.4.3  Celerity  and  Wave  Length 

In  order  to  further  understand  the  properties  of  internal  waves  at  the  lutocline  in  the 
Jiaojiang,  the  celerity  and  wave  length  are  calculated  here.  For  simplification,  the  flow  is 
treated  as  two-layered  system  with  fluid  densities  of  p,  and  P2  in  upper  and  lower  layers, 

respectively.  Once  the  reduced  gravity  g'  [ =g(p2-p,)/p]  and  the  lutocline  elevation  above 


the  bottom  are  known,  the  celerity  and  the  wave  length  can  be  calculated  according 


to  (Lamb,  1945) 


^tanh(A:CJ  , 


to: 


(4.22) 


where  k is  the  wave  number  ( =2-k/XJ.  From  the  measured  SSC  profiles  (Figure  4.9), 
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lutocline  elevation  (Figure  4.12)  and  the  modal  frequency  of  internal  waves  (Section  4.3.4), 
the  calculated  results  for  examples  (a),  (b)  and  (c)  in  Figure  4.12  are  presented  in  Table  4.3. 

For  these  calculations,  the  mean  lutocline  elevation  was  taken  for  each  ASSM  segment,  and 
an  iteration  method  was  used  in  solving  Eq.  (4.22). 

From  Eq.  (4.18),  one  can  calculate  the  critical  wave  length,  below  which  the 


interface  will  become  unstable,  i.e.,  the  Kelvin-Helmholtz  (Delisi  and  Corcos,  1973)  or 
Holmboe  (Browand  and  Wang,  1972)  instabilities  due  to  interfacial  shear.  By  equating  the 
two  sides  of  Eq.  (4.18)  one  obtains 


X = 


g(pl~p]) 


RL 


(4.23) 


Values  of  calculated  from  Eq.  (4.23)  are  given  in  Table  4.3.  Also  listed  in  Table  4.3  are 


the  Brunt- Vaisala  frequencies  calculated  from  Eq.  (4.7). 

From  Table  4.3  it  is  seen  that  the  high  frequency  waves  were  characteristically  in 
deep  water,  with  the  ratio  on  the  order  of  5 (»0.5),  whereas  the  low  frequency  waves 

were  close  to  the  shallow  water  regime,  with  on  the  order  of  0.07  (=0.05).  It  is  also 

observed  that  the  celerity  and  wave  length  increased  with  increasing  Richardson  number  for 
both  high  and  low  frequency  waves.  Observe  that  X^^  decreased  with  increasing  Richardson 

number  (ranging  from  2.71  m to  0.06  m).  Thus,  the  wave  lengths  for  high  frequency  internal 
waves  (ranging  from  0.39  m to  0.65  m)  are  between  the  maximum  and  minimum  critical 
wave  lengths  for  stability.  This  suggests  that  the  high  frequency  waves  are  generated  by 
forcing  due  to  interfacial  shear. 
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CHAPTER  5 

TURBULENCE  DAMPING  IN  FLUID  MUD 
5.1  Introduction 

As  stated  in  Chapters  1 and  2,  the  vertical  mixing  pattern  of  suspended  sediment  over 
the  lutocline  is  highly  dependent  on  the  nature  and  intensity  of  turbulence  in  both  layers. 
Laboratory  experiment  results  (e.g.,  Wolanski,  et  al.,  1989;  Mehta  and  Srinivas  1993; 
Winterwerp  and  Kranenburg,  1997)  and  field  investigations  (e.g.,  Jiang  and  Wolanski,  1998) 
have  revealed  that  upward  mixing  caused  by  the  instability  and  breaking  of  internal  waves 
at  the  lutocline  occurs  concurrently  Avith  turbulence  damping  within  the  fluid  mud  layer 
below. 

Turbulence  damping  by  suspended  sediment  was  early  examined  by  Einstein  and 
Chien  (1955).  They  argued  that  since  part  of  the  turbulent  energy  is  used  to  maintain  the 
sediment  particles  in  suspension,  turbulence  is  damped  by  the  suspended  sediment.  Since 
then,  turbulence  damping  in  the  fluid  mud  layer  and  consequent  (vertically)  asymmetric 
mixing  over  the  lutocline  have  been  commonly  reported  (e.g.,  Wolanski,  et  al.,  1992; 
Scarlatos  and  Mehta,  1993;  Kranenburg  and  Winterwerp,  1997;  Jiang  and  Wolanski,  1998). 
However,  there  remains  a lack  of  a theoretical  basis  as  well  as  any  direct  evidence  of  this 
feature  of  mixing  because  of  the  difficulty  in  observing  it  in  the  field. 
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Here,  turbulence  damping  in  the  fluid  mud  layer  is  examined  on  a phenomenological 
basis.  Its  effect  on  lutocline  formation  in  the  Jiaojiang  is  then  explained. 

5.2  Turbulence  Damning  and  its  Effect  on  Lutocline  Formation 
For  a simplified  treatment,  we  will  consider  a steady  uniform  flow,  treat  fluid  mud 
as  a Bingham  plastic  and  omit  advective  effects. 

The  rate  of  production  of  turbulent  energy  associated  with  the  Reynolds  stress  (per 
unit  volume)  can  be  expressed  as  x^^du/dz  (Rossby  and  Montgomery,  1935).  Based  on  the 

Prandtl  mixing  theory,  Einstein  and  Chien  (1955)  proposed  following  formula  for  the 
Reynolds  stress,  t^,  in  a uniform  flow 


— m'(P,-c)w-m'c(w-o)  ) 


P.Po  , 


Po^- 


du 
dz , 


(5.1) 


where  u'  is  the  turbulent  fluctuation  of  the  horizontal  velocity  in  the  x direction,  / is  the 

momentum  mixing  length,  and  the  overbars  denote  time-averaging.  Hence  the  rate  of 
production  of  turbulent  energy  becomes 


du 
' dz 


P,Po 


^41)' 


(5.2) 


The  rate  of  work  done  against  buoyancy  due  to  stratification  is  (Odd  and  Rodger, 


1978) 


dz 


du 

du 

dz 

fn  fn 

dz 

(5.3) 
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where  is  the  vertical  flux  of  buoyancy  and  is  the  ratio  of  the  mass  mixing  length,  /^, 

to  the  momentum  mixing  length,  /^,  i.e.,  the  turbulent  Schmidt  number. 

The  rate  of  work  done  against  gravity  is  (Hunt,  1954) 


(5.4) 

where  A/  is  the  vertical  mass  flux  and  c'is  the  turbulent  fluctuation  of  SSC.  For  a steady 

flow,  from  the  sediment  and  water  continuity  considerations  and  the  Prandtl  mixing  theory, 
Einstein  and  Chien  (1955)  obtained 


w'c'=- 


p -c 


-CO). 


(5.5) 


Substituting  Eqs.  (3.27)  and  (5.5)  into  Eq.  (5.4),  rate  of  work  done  against  gravity  becomes 


g(p,-Po)(p,-c) 

gM  ^ : cw. 


2 

p. 


(5.6) 


The  rate  of  work  done  against  cohesion  and  interactions  between  the  floes  in  the  fluid 


mud  is 


^ du  _/  \du 


dz 


(5.7) 


where  is  the  total  shear  stress  due  to  cohesion  and  interactions  between  the  floes  and 

is  the  shear  stress  due  to  the  interactions  between  floes.  Bagnold  (1954;  1956)  examined  the 
normal  and  tangential  stresses  in  granular  flows  and  suggested  that  they  may  be  expressed 
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as  functions  of  the  shear  rate,  du/dz,  and  a “linear  sediment  concentration”, 
characterizes  the  relative  surface  proximity  between  sediment  particles  (Figure 
related  to  the  sediment  concentration,  c,  by 

1 

where  is  the  maximum  concentration  corresponding  to  grain-grain  contact. 


c^,  which 
5.1).  is 

(5.8) 


< > 


Figure  5.1.  Definition  of  linear  sediment  concentration,  c^,  and  its 
relationship  with  sediment  concentration,  c.  the  floe  diameter. 


Observe  that  increases  drastically  as  c approaches  (Figure  5.1).  For  small. 


light  grains  in  a very  viscous  fluid,  Bagnold  (1954;  1956)  found  that  the  behavior  of  the 
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mixture  of  fluid  and  cohesionless  grains  is  dominated  by  viscosity  due  to  interactions 
between  particles,  and  termed  it  the  macro-viscous  regime.  In  this  regime  the  shear  stress  has 
the  form 


tancj) 


V = 1.3(l+c,) 


1+-A 


/ 


du 


(5.9) 


where  (})^  is  the  dynamic  angle  of  repose,  which  was  found  to  be  (J)^=37°  (Bagnold,  1954; 

1956),  and  p is  the  dynamic  viscosity  of  the  fluid.  To  extend  the  applicability  of  this  formula 
to  fluid  mud,  it  is  reasonable  to  assume  that  the  behavior  of  a mixture  of  fluid  and  cohesive 
sediment  is  phenomenologically  analogous  to  the  behavior  of  a mixture  of  fluid  and 
cohesionless  grains.  Accordingly,  Eq.  (5.9)  can  be  rendered  applicable  to  the  fluid  mud 
merely  by  taking  the  dynamic  viscosity  p to  be  that  of  the  fluid  mud,  and  treating  floes  as 
basic  particles. 

Following  the  original  arguments  of  Rossby  and  Montgomery  (1935),  under  an 
assumed  equilibrium  condition  the  sum  of  the  kinetic,  potential  and  dissipated  energy  in  a 
stratified  and  cohesive  flow  per  unit  mass  can  be  considered  to  be  the  same  as  that  in  a 
homogeneous,  non-cohesive  flow  at  identical  shear  rates.  Thus,  combining  Eqs.  (5.2),  (5.3), 
(5.5),  (5.7)  and  (5.9)  one  obtains 


2 


dz  I 


P.Po 


Po^i 


li  dwV 

ATz 


Sr  mm  ^ ^ 

dz  dz 


. ^(P.-PoXP.-^^)  , X ^ 

+ C0)^+I.3tan4)/1  +q) 

P. 


f 

1+^ 

1 2j 

du 
^ dz 


(5.10) 
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Solving  Eq.  (5.10)  for  /pleads  to 


(5.11) 


with 


P.Po 

; _ (P.~P0)(P,-g)^ 

pjpo 


i?/^  = 1.3tanc|)j(l  +q)(l  +-^) ti 


(5.12) 


Ri^ 


9(floidu/dzf 


Here  is  a dimensionless  variable  dependent  on  SSC,  Ri^  is  the  ratio  of  the  potential 
energy  of  the  sediment  settling  flux  to  the  production  of  turbulent  energy,  ll^{du/dzf , 

(denoted  as  F^),  Ri^  is  the  ratio  of  the  viscous  force  due  to  the  interactions  between  floes  in 
the  fluid  mud  to  the  Reynolds  stress,  and  Rig  is  the  ratio  of  Bingham  yield  stress  to  the 
Reynolds  stress. 

Equation  (5.1 1)  shows  that  at  equilibrium  the  higher  the  potential  settling  flux  of 
suspended  sediment  (as  would  occur  in  quiescent  water),  the  smaller  the  turbulent  mixing 
length.  In  other  words,  the  turbulent  kinetic  energy  is  damped  because  part  of  it  is  used  to 
maintain  floes  in  suspension.  Due  to  hindered  settling,  the  sediment  settling  flux  will 
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decrease  with  increasing  SSC  (Ross  and  Mehta  1989).  Thereafter,  cohesion  and  the 
interactions  between  floes  can  be  expected  to  play  an  increasingly  important  role  in 
turbulence  damping.  To  demonstrate  this  behavior,  one  may  introduce  the  following 
formulas  and  parameters: 

(1)  Settling  Velocity:  Here  the  floe  settling  velocity  formula  (3.30)  is  applied,  with 
a=0.085,  6=10  kg  m'^  a=1.5,  P=1.6,  and  without  considering  the  effects  of  salinity  and 
temperature  on  flocculation. 

(2)  Bingham  Yield  Strength:  The  Bingham  yield  strength  is  expressed  as  a function 
of  SSC  as  follows 

(5.13) 

Here  and  are  sediment-dependent  coefficients.  Following  Owen  (1970)  and  Odd  and 
Rodger  (1986),  ag=7.36xl0“^  and  P^=2.33  are  chosen. 

(3)  Momentum  Mixing  Length:  In  the  near-bottom  layer  of  a homogeneous,  non- 
cohesive  flow,  the  following  Prandtl  (1925)  approximation  for  the  momentum  mixing  length 
can  be  applied 

(5.14) 

where  is  the  elevation  above  the  bottom,  and  von  Karman  constant  k is  normally  taken 
as  0.4. 

(4)  Velocity  Profile:  In  the  near-bottom  layer,  the  velocity  distribution  will  be 
considered  to  be  a locally  logarithmic,  i.e.. 
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“.n 

w(z.)=— In— 


K 


(5.15) 


Where  is  the  bottom  frictional  velocity  in  homogenous  flow. 

(5)  Viscosity  of  Fluid  Mud:  Following  Odd  and  Rodger  (1986),  the  viscosity  of  fluid 
mud,  |i,  with  SSC<80  kg  m ^ , which  is  the  critical  SSC  for  soil  formation,  will  be  taken  as 
0.01  Pa  s. 

(6)  Water  and  Granular  Densities:  Values  of  water  and  granular  densities  taken  here 
are  Pq=1,000  kg  m'^  and  =2,650  kg  m■^ 


(7)  Maximum  Concentration:  As  mentioned  above,  for  interactions  between  particles 
in  fluid  mud  the  floes  are  treated  as  basic  particle  units.  Hence  the  maximum  concentration, 
^max’  should  Correspond  to  floc-floc  contact.  Floes  tend  to  be  very  loosely  bound  and  are 

light  in  water,  with  a typical  bulk  density,  p^,  of  about  1,080-1,150  kg  m which  will 

contain  nearly  95%  by  volume  of  water  locked  within  the  interstitial  particulate  fabric  (Mehta 
and  Li,  1997).  Accordingly,  the  maximum  concentration  corresponding  to  the  floc-floc 
contact  can  be  expressed  as 


= 


(Pp-P/K 

C„ 


(5.16) 


Where  is  the  maximum  volumetric  floe  content  corresponding  to  floc-floc  contact.  In 


analogy  with  the  water-sand  mixture  (Bagnold,  1956),  c^^=0.65  will  be  chosen. 
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Substituting  the  above  parameters  and  expressions  into  Eq.  (5.1 1),  the  ratio  of  / 

is  calculated  and  plotted  in  Figure  5.2  as  a function  of  SSC  and  the  production  of  turbulent 
energy,  F^,  below  SSC<80kg  m ^ without  considering  the  effects  of  buoyancy.  Also  shown 

is  the  sediment  settling  flux  as  a function  of  SSC.  It  is  observed  that  within  the  flocculation 
settling  range,  maximum  turbulence  damping  is  consistent  with  the  maximum  settling  flux. 
Thus,  a lutocline  can  be  expected  at  the  elevation  of  the  maximum  settling  flux  in  the  water 
column.  As  soon  as  the  lutocline  is  formed,  the  buoyancy  effect  due  to  sediment-induced 
stratification  will  enhance  turbulence  damping  as  indicated  in  Eq.  (5.1 1).  This  in  turn  will 
accentuate  the  lutocline.  Figure  5.2  also  shows  that  in  the  range  of  hindered  settling, 
turbulence  damping  decreases  due  to  decreasing  sediment  settling  flux  for  F <0.0001 

m^  s However,  above  SSC>40-50  kg  m damping  increases  drastically  as  the  effects 

of  cohesion  and  interactions  between  floes  become  significant.  Turbulence  collapse 
subsequently  occurs.  Thus,  turbulence  damping  is  governed  by  the  potential  settling  flux  in 
the  flocculation  settling  range,  and  by  cohesion  and  interactions  between  floes  in  the 
hindered  settling  range  with  SS040-50  kg  m Overall,  flocculation  settling,  turbulence 

damping  due  to  settling  flux,  cohesion  and  interactions  between  floes  in  fluid  mud  as  well 
as  the  buoyancy  effect  due  to  sediment-induced  stratification  contribute  to  the  formation  of 
a lutocline.  This  analysis  thus  replicates  the  conclusion  arrived  at  by  Ross  and  Mehta  (1989) 
based  solely  on  the  concept  of  density  stratification.  However,  the  explanation  provided  here 
introduces  new  parameters,  namely  Ri^,  Ri^  and  Ri^. 
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Settling  flux,  Cl»,cx10^  (kg  s ') 

0.002  0.004  0.006  0.008  0.01 


Figure  5.2.  Relative  momentum  mixing  length  calculated  from  theoretical  formula  (solid 
lines)  and  field  data  (data  points)  and  settling  flux  (dashed  line)  as  fimctions  of  SSC. 


Also  observed  in  Figure  5.2  is  that  turbulence  damping  is  dependent  not  only  on  SSC, 
but  also  on  flow.  Thus,  beyond  F^>~0.0002  m ^ s the  ratio  / approaches  unity,  i.e., 
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there  is  no  significant  damping  and,  as  a result,  mixing  will  occur  in  two  ways,  i.e.,  both 
upward  and  downward.  The  lutocline  will  be  ruptured  and  become  less  distinct.  This 
conclusion  is  supported  by  the  observations  of  Jiang  and  Wolanski  (1998)  in  the  Jiaojiang, 
where  it  was  foimd  that  a distinct  lutocline  only  appears  during  slack  water  at  neap  tide  when 
the  production  of  turbulent  energy,  F^,  is  usually  less  than  0.0002  m^  s . 


To  examine  the  relationship  between  lutocline  formation  and  the  flow  condition,  an 
analysis  of  vertical  profiles  of  SSC  and  velocity  in  the  1991  and  1994  field  campaigns  in  the 
Jiaojiang  (Table  4.1)  is  considered  here.  The  result  is  shown  in  Figure  5.3  in  terms  of  a 
lutocline  strength  index,  as  a function  of  the  production  of  turbulent  energy.  Here,  is 

defined  as 


[dc/  dz) 

r _ ' ' 'mean 

where  subscripts  “max"  and  “mean"  denote  the  maximum  and  the  mean  values  of  the  SSC 
gradient,  dc/ dz , over  the  water  column,  respectively.  = 1 indicates  that  there  is  a uniform 


vertical  distribution  of  the  SSC  and  no  lutocline.  The  vertical  gradient  of  SSC  is  calculated 
from 


dc  _^k*\ 
dz  dz 


(5.18) 


where  subscript  “Z”  denotes  the  measurement  elevation  and  dz  is  the  height  between  k and 
Z+l.  Figure  5.3  shows  the  results  and  the  mean  trend  line.  It  is  evident  that  the  lutocline 
becomes  less  distinct  as  the  turbulent  energy  production  F^>~ 0.0003  m^  s Relatively 


125 


high  value  of  F^,  at  which  the  lutocline  became  less  distinct,  are  found  because  the  buoyancy 
effect  due  to  sediment-induced  stratification  was  not  considered  in  Figure  5.2. 


Production  of  turbulent  energy,  F,  (m^  s'^) 


Figure  5.3.  Lutocline  strength  index  as  a fimction  of  turbulence 
energy  production  based  on  measured  profiles  of  SSC  and 
velocity.  The  equation  represents  the  best-fit  line. 


5.3  Mixing  Length  in  the  Jiaoiiang 

From  the  Jiaojiang  ASSM  field  data  (e.g.,  Figure  4.10),  it  is  possible  to  examine  the 
standard  deviation  of  SSC,  o^,  which  has  the  form 
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^4 


'^d(  \ 2 


^-1 
V c ) 


dt 


(5.19) 


where  T^  is  the  period  over  which  is  calculated  and  the  overbar  expresses  time  averaging 

over  a period  of  T^.  Here  7)^  will  be  taken  as  2.5  min  as  a characteristic  duration,  c is 

calculated  from  Eq.  (4.1)  by  using  ASSM  records.  Considering  that  Eq.  (5.19)  is  applied  to 
a fixed  elevation  above  bottom,  i.e.,  /i'=constant,  Eq.  (4.1)  becomes 


c=AF„ 


(5.20) 


where  A is  a sediment  and  elevation  dependent  constant.  By  assuming  that  is  proportional 

to  the  local  Prandtl  mixing  length  and  the  vertical  gradient  of  SSC,  one  obtains 

( -\2 


dc 
\ dzj 


(5.21) 


For  a homogeneous  flow,  assuming  that  the  shear  stress  near  the  bottom,  x , is  constant,  the 


local  Prandtl  mixing  length  is  obtained  as 


Co 


\ {du/dzf  du/dz 
Combining  Eqs.  (5.19),  (5.20),  (5.21)  and  (5.22),  one  obtains 


(5.22) 


du  

C dz  P'a 


Co  Pm".  BFl/dz\ 


±r^-i 

T J ^2 


dt 


(5.23) 


127 


UCS  of  ^ rrJ^ mO  were  evaluated  by  Eq.  (5.23)  using  ASSM  records  and  measured 

velocity  profiles  (Figure  4.6).  The  results  are  plotted  in  Figure  5.2,  where  it  is  observed  that 
the  trends  in  mixing  length  follow  Eq.  (5.1 1). 

5.4  Modified  Vertical  Momentum  and  Mass  Diffusion  Coefficients 


The  above  results  suggest  that  to  calculate  the  turbulent  diffusion  coefficients  in  a 
cohesive  sediment  transport  model,  it  is  essential  to  consider  the  effects  of  potential  settling 
flux,  cohesion  and  interactions  between  floes  as  well  as  sediment-induced  stratification,  as 
embodied  in  Eq.  (5.1 1).  Accordingly,  the  turbulent  diffusion  coefficient  formulas  (3.39)  and 
(3.40)  are  modified  following  Munk  and  Anderson  (1948)  analog  as: 

Momentum  diffusion  coefficient: 


^v=^v0 


1 


a.+b,Ri 

fn  \ 


+A 


vb 


(5.24) 


Mass  diffusion  coefficient: 


a+b.Ri 

tn  I 


+K. 


vb 


(5.25) 


where  (<1)  and  (<1)  are  the  sediment-dependent  coefficients  relating  the  effect  of 


suspended  sediment  on  the  mixing  length. 


CHAPTER  6 

LUTOCLINE  DYNAMICS  IN  THE  JIAOJING 
6.1  Introduction 

COHYD-UF  and  COSED-UF  are  now  applied  to  simulate  lutocline  and  associated 
fluid  mud  dynamics  in  the  Jiaojiang.  The  modeled  results  include  tidal  variations  of  velocity, 
SSC  and  the  lutocline  layer,  lutocline  layer  thickness  and  hysteresis  loops  of  SSC,  and  these 
are  compared  with  the  data  obtained  in  situ.  Relevant  formulations  of  the  flow  and 
sedimentary  processes  modeled  are  summarized  in  Section  6.2.  Applications  of  the  numerical 
models  are  described  in  Section  6.3.  Finally,  simulations  along  with  the  data  are  presented 
in  Section  6.4. 

6.2  Parameters  for  Flow  and  Sedimentary  Processes 
Descriptions  of  the  models  and  relevant  flow  and  sedimentary  formulations  are  given 
in  Chapters  3, 4 and  5.  Table  6.1  summarizes  these  formulations  and  relevant  parameters  in 
modeling  sediment  dynamics  in  the  Jiaojiang.  The  most  important  parameters,  including  the 
settling  velocity,  erosion  rate  constant,  momentum  and  mass  diffusion  coefficients, 
consolidation  rate,  and  the  effective  roughness  of  the  bed,  were  determined  through  analyses 
of  field  data  and/ or  model  calibration.  Others  were  introduced  from  the  works  of  previous 
researchers. 
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Table  6.1.  Flow  and  sedimentary  process  formulations  and  parameters 
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Table  6.1— continued 
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Table  6.1 --continued 
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6.3  Model  Application 
6.3.1  Modeled  Domain.  Initial  and  BoundatY  Conditions 

The  modeled  region  includes  the  area  from  the  mouth  to  the  upstream  tributary 
(Figure  4.1),  with  a length  of  13.4  km.  The  domain  was  discretized  using  spatial  steps 
m and  ao=0.1,  with  the  total  number  of  rectangular  cells 

=67x13x10=8,710.  Figure  6.1  shows  the  bathymetry  and  numerical  mesh  in  the  horizontal 
plane  within  the  modeled  domain.  A time  step  a/=30  s was  adopted  based  on  stability 
constraints  resulting  from  the  numerical  scheme  involving  the  back-tracing  approach 
mentioned  in  Section  3.5.4  [Eq.  (3.50)].  The  lowest  layer  near  the  bottom  was  0.05// above 
the  bed  and  the  highest  layer  near  the  surface  was  0.05//  below  the  instantaneous  water 
surface. 

The  initial  conditions  for  COHYD-UF  runs  were  as  follows:  (a)  All  velocities  were 
set  equal  to  zero,  and  (b)  water  elevation  at  each  point  was  generated  through  linear 
interpolation  of  the  prescribed  values  at  the  open  boundaries  based  on  measurements.  For 
COSED-UF,  the  initial  SSC  was  obtained  by  the  same  interpolation  approach  using  the 
values  at  the  open  boundaries  from  measurements.  Once  a stable  flow  field  resulted  from 
COHYD-UF  model,  which  took  about  6 hr  (one-half  tidal  cycle),  COSED-UF  run  was 
initiated. 

Stations  T1  and  T5  (Figure  4.1)  were  selected  to  prescribe  the  open  boundary 

conditions  for  the  water  surface  elevation.  The  following  1 1 tidal  constituents  were 

considered  at  these  boundaries:  Q,,  Oj,  P,,  Kj,  N2,  M2,  S^,  K^,  M^,  M^  and  M . All 

6 ^6 
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harmonic  components  were  assumed  to  be  constant  in  phase  and  amplitude  over  the  cross- 
section  of  each  open  boundary. 


(a) 


Upstream 


dx=dy=2QQ  m 
— ►x 


Mouth 


(b) 


Figure  6.1.  Bathymetry  in  the  modeled  domain  of  the  Jiaojiang  (a),  where 
the  datum  is  mean  water  level  and  the  regions  enclosed  within  dotted  lines 
are  mudflats,  and  the  numerical  plane  mesh  in  the  horizontal  plane  (b). 


Measured  SSC  values  at  sites  Cl  and  C3  (Figure  4.1)  were  inputted  as  the  open 
boundary  conditions  for  COSED-UF.  Since  there  were  only  one  measurement  station  over 
the  cross-section  of  each  open  boundary  (located  in  the  middle  as  shown  in  Figure  4.1),  it 
was  assumed  that  the  SSC  was  uniform  over  the  entire  cross-section  of  the  open  boundaries. 
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In  reality  there  is  always  a certain  amount  of  lateral  variation  of  SSC,  even  in  narrow 
channels.  Consequently,  this  uniformity  approximation  can  lead  to  a noticeable  error  in 
simulations  results  within  the  modeled  domain.  This  aspect  will  be  discussed  further  in 
Section  6.4.  Because  SSC  observations  were  made  at  variable  elevations  using  a turbidimeter 
lowered  from  the  boat,  the  SSC  values  at  fixed  elevations  had  to  be  obtained  through  cubic 
spline  interpolations. 

Currents  and  SSC  at  sites  C2  and  C4  (Figure  4.1)  were  used  to  verify  the  simulations. 
Model  tests  were  run  over  the  same  period  as  the  field  campaign,  i.e.,  from  November  4, 
1994  (spring  tide)  to  November  12, 1994  (neap  tide).  Simulated  flow  field  and  SSC  profiles 
at  the  same  time  as  observations  at  C2  and  C4  were  outputted  for  verification  and  analysis. 
6-3.2  Sediment  Deposition.  Erosion.  Consolidation  and  Entrainment 

1 . Deposition  and  Erosion:  When  sediment  deposition  occurs  (Section  3.3.3 

and  Table  6.1).  Within  a time  step,  At,  the  rate  of  deposition,  m^,  is  assumed  to  be  constant 

and  calculated  by  the  formula  in  Table  6.1.  When  where  is  a function  of  sediment 

concentration,  bottom  erosion  occurs  (Section  3.3.3  and  Table  6.1).  First  eroded  is  the 
consolidating  layer  and  then  the  fully  consolidated  bed.  Within  time  step.  At,  erosion 
continues  until  is  encountered.  When  it  is  assumed  that  deposition  and 

erosion  are  in  equilibrium,  i.e.,  the  sediment  vertical  flux  at  the  bed-fluid  interface  m -m , 

is  zero.  Deposition  and  erosion  are  incorporated  in  COSED-UF  through  the  bottom  boundary 
condition,  i.e.,  Eq.  (3.17). 
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2.  Consolidation:  The  deposited  material  is  combined  with  the  consolidating  sediment 
layer  having  an  initial  sediment  concentration  of  (Table  6.1),  and  goes  through 

consolidation  (Sections  3.3.3  and  3.5.3).  When  sediment  concentration  near  the  bottom  of 
the  consolidating  layer  is  greater  than  the  maximum  compaction  concentration,  this 

portion  of  the  deposit  is  combined  with  the  fully  consolidated  bed. 

During  deposition,  bottom  erosion  and  consolidation  of  freshly  deposited  sediment, 
although  the  heights  of  the  consolidating  and  the  fully  consolidated  layers  vary  with  time,  the 
bottom  datum  remains  constant. 

3.  Entrainment:  Interfacial  entrainment  discussed  in  Section  3.3.3  and  turbulent 
diffusion  in  Section  3.4.2  and  Chapter  5 are  facets  of  interfacial  mixing  in  a stratified  flow, 
and  are  referred  to  as  pure  entrainment  mixing  and  the  pure  turbulent-diffusion  mixing. 
Grubert  (1990)  noted  that  pure  entrainment  mixing  takes  place  when  the  interfacial  transition 
layer  (or  lutocline  layer)  is  in  a subcritical  state.  In  this  condition,  cusp  motions  generated  by 
interfacial  instabilities,  transfer  volumes  of  fluid  in  either  direction  to  a greater  or  lesser 
extent,  depending  on  turbulence  in  each  layer. 

Pure  turbulent-diffusion  mixing  occurs  when  the  transition  layer  is  in  a supercritical 
state.  In  this  condition,  violent  vortex  motions  exchange  equal  volumes  of  fluid  between 
layers.  In  reality,  these  two  mixing  processes  are  usually  superimposed  on  each  other,  and 
the  critical  condition  is  dependent  on  local  site-specific  processes.  By  reexamining  the 
experimental  data  of  previous  researchers  including  Ellison  and  Turner  (1959),  Lofquist 
(1960),  Kato  and  Phillips  (1969),  Moore  and  Long  (1971),  Wu  (1973),  Chu  and  Vanvari 
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(1976),  Kantha,  et  al.,  (1977),  Bo  Pedersen  (1980),  Kit,  et  al.,  (1980),  Buch  (1981),  and 
Narimousa,  et  al.,  (1986),  Christodoulou  (1986)  found  that  in  the  Richardson  number  range 
0.01  </?/q:£  1 the  interface  is  turbulent  and  mixing  is  produced  directly  by  eddies  in  the  mixed 

layer.  In  the  range  0.1  10,  mixing  is  less  turbulent  and  is  produced  by  Kelvin- 

Helmholtz  instabilities.  In  the  range  1 100,  where  shear  is  weak,  mixing  is  produced 

by  intermittent,  cusp-generated  entrainment.  Finally,  when  /?/q^100,  mixing  is  due  to 

molecular  diffusion.  Mehta  and  Srinivas  (1993)  conducted  experiments  of  fluid  mud 
entrainment  in  a shear  flow  using  a race-track  flume.  They  observed  that  for  Ri^<5, 

interfacial  entrainment  appeared  to  be  turbulence  dominated.  As  Ri^  increased  above  ~5,  the 

interface  become  convoluted  with  large,  irregular  undulations.  Entrainment  appeared  to  be 
dominated  by  interfacial  wave  breaking  in  which  wisps  of  fluid  were  episodically  ejected  into 
the  mixed  layer. 

Thus,  considering  the  fact  that  at  very  low  Richardson  number,  i.e.,  Ri^->0,  mixing 

over  the  lutocline  is  dominated  by  turbulent  diffusion,  and  in  this  situation  the  entrainment 
formula  (Table  6.1)  becomes  inapplicable  due  to  The  model  only  activates  the 

entrainment  formula  when  Ri^^  1 . 

6.4  Flow  and  Sediment  Dynamics 

6.4.1  Flow  Field 

Spring  tide  peak  flows  and  slack  water  within  the  modeled  domain  are  shown  in 
Figures  6.2-6.5.  The  flow  field  shows  the  following  noteworthy  features: 
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(1)  During  peak  flow  the  velocity  decreases  gradually  from  the  mouth  to  upstream, 
while  during  slack  water  it  decreases  from  upstream  to  the  mouth. 

(2)  Maximum  velocities  occur  along  the  channel  center,  in  the  channel  segment  near 
convex  shorelines,  i.e.,  shorelines  protruding  into  the  water  area,  and  at  other  narrow 
sections. 

(3)  Near  concave  shorelines  the  velocities  are  lower.  There  mudflats  are  usually  found 
[Figure  6.1(a)]. 

6.4.2  Tidal  Variation  of  Velocity 

Model  simulated  tidal  velocities  are  shown  in  Figures  6.6(a)-6.9(a).  The  flood 
current  was  stronger  than  the  ebb  current.  For  example,  at  0.4// elevation  the  flood  peak  at 
C2  during  spring  tide  was  around  1.70  m s , and  the  ebb  peak  around  1.08  m s''  [Figure 

6.6(a)],  with  a ratio  of  1.57.  However,  the  peak  ebb  flow  lasted  considerably  longer  than  the 
peak  flood  flow.  From  Figures  6.6(a)-6.9(a)  it  is  found  that  the  duration  of  peak  ebb  was 
around  5.34  hr,  whereas  that  of  the  peak  flood  was  around  2.05  hr,  with  a ratio  of  2.60. 

A noteworthy  difference  between  simulations  and  observations  is  that  the  model 
under-predicted  the  ebb  peak  at  C4  [Figures  6.8(a)  and  6.9(a)].  To  highlight  this  difference 
further,  vertically-averaged  maximum  velocities  at  C2  and  C4  are  given  in  Table  6.2.  It  is 
evident  that  the  greatest  difference  between  simulations  and  observations  occurred  at  C4 
during  peak  ebb,  with  a difference  of  about  0.37  m s . At  C4  the  measured  peak  ebb  was 

greater  than  peak  flood  in  contrast  to  other  sites  where  flood  was  dominant.  Table  6.2  also 
shows  that  the  flood  current  at  C2  was  greater  than  at  C4,  whereas,  ebb  at  C4  was  greater 
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Figure  6.6.  Tidal  velocity  at  0.4//  (a),  and  SSC  at  0.25//  (b)  Figure  6.7.  Tidal  velocity  at  0.4//  (a),  and  SSC  at  0.25//  (b) 

and  at  0.75//  (c)  at  C2  during  a spring  tide.  Solid  lines  are  and  0.75//  (c)  at  C2  during  a neap  tide.  Solid  lines  are 

simulations  and  dashed  lines  represent  field  data  collected  simulations  and  dashed  lines  represent  field  data  collected 

during  1800  hr,  1 1/05/94  to  1700  hr,  1 1/06/94.  during  1000  hr,  1 1/10/94  to  0900  hr,  1 1/1 1/94. 


141 


o 


~a  ^ 

§ 5 

Vi  O 
- C ^ 

. o 


52  ^ 

in  3 

d.§ 

^ </3 

c3  p 

U S 

cn  “ 

T3  S 
§ T3 


3 


cn 


60 

c 

•c 

3 

T3 

•o 

(U 

4-^ 

o 

ID 


ID 

Tt  T3 

d 


C6 


O 

_o 

13 

> 

13 

T3 


(D 

C 

c6 

60 

G 

■c 

G 

T3 

u 


c3 


O 


o\ 

\D 

G a; 

60  in 
£ 

o 


•a 

u 

J3 


■It 


H;  (u  = 


o 

o 

cd 

cd 

•O 

2 

13 

<4-l 

c 

o 

c/i 

H 

a, 

<u 

c 


u 

x: 


ON 

o 


rsi 


i/n 


1 

1 

\ ^ 

\ 

\ 

x\ 

1 1 



/ 

\ w 
! 

^ y\ 

\ 

A....1... 

'n\ 

1 i / 

to 

C'J 


g 


wn 

CS 


cd 

U 

C/D 

C/3 

•o 

§ 


§ 

c« 

ID 

C 


T3 

U 

4-* 

o 

o ^ 
o ^ 
G iB 
■“  P 


S c6 
O T3 
C/3  -o 

oJ  13 

T3  C 


60 

C 


« -a  ~ 


a; 

d 

'S 

4— > 

‘o 

_o 

13 

> 

•o 


l-l 

CD 

Vi 

G 


SP 

c 


T}- 

pi 


G 

T3 

u 

t3 


u m 
•S  5 

T3  £ 

x:  uT 


"O 

T3 

§ 


o s 2 


(,.«  in)  iC}|Ooi0A  (^.IH  Sji)  OSS 


(e.™  3^1)  OSS 


00 

d 

a 

3 

60 

E 


a; 

m 

r~ 

d 

T3 

§ 


60 

_C 

'C 

3 

XJ 


142 


than  at  C2.  These  observations  colleetively  suggest  that  there  was  a horizontal  eirculation  in 
the  clockwise  sense  wdthin  the  domain.  This  effect  is  discussed  further  in  Section  6.4.4. 


Table  6.2.  Vertically-averaged  maximum  velocities  at  sites  C2  and  C4  (Unit:  ms') 


Site 

Tide  range 

Observation 

Simulation 

Flood 

Ebb 

Flood 

Ebb 

r? 

Spring  tide 

1.60 

1.42 

1.66 

1.40 

Neap  tide 

1.33 

0.99 

1.56 

1.16 

C4 

Spring  tide 

1.24 

1.38 

1.26 

1.03 

Neap  tide 

1.16 

1.21 

1.36 

0.82 

6.4.3  Tidal  Variation  of  SSC 

The  tidal  variation  of  SSC  near  the  bottom  (specifically  0.25// above  the  bottom)  and 
near  the  surface  (0.25// below  the  surface)  are  presented  in  Figures  6.6(b)-6.9(b)  and  6.6(c)- 
6.9(c),  respectively.  The  following  noteworthy  differences  in  the  magnitude  of  SSC  between 
simulations  and  observations  occur: 

(1)  Near  the  bottom,  the  model  mostly  under-predicts  SSC  during  flood  and  over- 
predicts it  during  ebb  [Figures  6.6(b)-6.9(b)]. 

(2)  Near  the  surface  the  model  mostly  under-predicts  SSC  over  the  entire  tidal  cycle 
[Figures  6.6(c)-6.9(c)]. 

Tidal  variations  of  SSC  show  the  following  characteristics: 

(1)  Variations  of  SSC  with  Tidal  Range:  SSC  exhibits  quantitatively  the  same 
behavioral  pattern  during  spring  and  neap  tides.  However,  near  the  bottom,  SSC  during 
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spring  tide  is  less  than  at  neap  [Figures  6.6(b)-6.9(b)],  being  usually  10-20  kg  m during 
spring  and  10-25  kg  m during  neap.  In  contrast,  near  the  surface,  SSC  during  spring  is 
greater  than  during  neap  [Figures  6.6(c)-6.9(c)],  being  usually  5-10  kg  m during  spring 
and  2-5  kg  m during  neap. 

(2)  Vertical  Variations  of  SSC:  As  might  be  expected,  SSC  near  the  bottom  was 
always  greater  than  near  the  surface.  Near  the  bottom  it  was  seldom  less  than  10  kg  m 

[Figures  6.6(b)-6.9(b)],  whereas,  near  the  surface  it  was  seldom  greater  than  10  kg  m 
[Figures  6.6(c)-6.9(c)]. 

(3)  Variations  of  SSC  during  Flood  and  Ebh:  During  high  water  slack,  “low”  SSC, 
which  is  defined  here  as  the  value  at  the  “trough”  of  SSC  tidal  variation,  occurred  both  near 
the  bottom  and  the  surface,  with  a lag  of  0.5-1  hr  near  the  surface.  During  low  water  slack 
peak  SSC,  which  is  defined  in  an  analogous  way  as  the  value  at  the  crest  of  SSC  tidal 
variation,  occurred  near  the  bottom  whereas  low  SSC  occurred  near  the  surface.  During  peak 
ebb,  peak  SSC  occurred  over  the  entire  water  column.  During  peak  flood,  peak  SSC  occurred 
near  the  surface  whereas  low  SSC  occurred  near  the  bottom. 

The  behavior  of  SSC  reflects  cumulative  effects  on  SSC  due  to  floe  settling,  tidal 
current  asymmetry  (with  a stronger  flood  peak  and  a weaker,  longer  duration  of  ebb  peak), 
entrainment  of  the  lutocline,  turbulence  damping  in  the  fluid  mud  layer  and  the  large  ratio 
of  tidal  range  to  mean  water  depth. 
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During  slack  water,  due  to  comparatively  low  vertical  turbulent  diffusion  over  the 
water  column,  higher  turbulence  damping  in  the  fluid  mud  layer  and  decreasing  mixing  over 
the  lutocline,  sediment  tends  to  settle  and  concentrate  near  the  bottom  where,  as  a result,  high 
SSC  occurs  in  the  hindered  settling  range.  As  a further  result,  peak  SSC  occurs  near  the 
bottom  and  low  SSC  near  the  surface,  along  with  a thick  fluid  mud  layer  and  a stable 
lutocline  (Odd  and  Rodger,  1986;  Wolanski,  et  al.,  1988).  This  was  the  situation  during  slack 
water  in  the  Jiaojiang  with  the  exception  of  the  condition  near  the  bottom  during  high  water 
slack.  It  is  believed  that  the  large  ratio  of  the  tidal  range  to  the  water  depth  causes  this  latter 
phenomenon.  As  stated  in  Section  4.1,  the  mean  tidal  range  in  the  studied  domain  was  about 
4 m and  the  mean  water  depth  (below  MSL)  about  4-7  m with  a ratio  of  about  0.8  between 
tidal  range  and  water  depth.  Thus,  during  flood  the  water  depth  increased  dramatically.  This 
increase  plus  advection  of  lower  SSC  water  from  the  region  beyond  the  modeled  domain 
caused  the  SSC  near  the  bottom  to  be  diluted.  As  shown  in  Figure  6.10,  the  vertically- 
averaged  SSC  at  Cl  during  flood  was  always  less  than  that  at  C2  and  C4,  with  a difference 
of  about  4-13  kg  m , which  accounts  for  about  50%-85%  of  the  vertically-averaged  values 

of  SSC  at  C2  and  C4.  It  is  this  difference  plus  the  large  ratio  of  tidal  range  to  water  depth  that 
causes  dilution  during  flood.  It  also  implies  that  sediment  transport  seems  to  be  influenced 
strongly  by  advection. 
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(a) 


(b) 


Figure  6.10.  Time  series  of  vertically-averaged  SSC  during  a spring  tide  (a)  and  a neap 
tide  (b).  During  spring  tide,  the  data  at  sites  Cl  and  C3  began  at  1700  hr,  1 1/04/94 
and  at  sites  C2  and  C4  at  1800  hr,  1 1/05/94.  During  neap  tide,  the  data  at  sites  Cl 
and  C3  began  at  1000  hr,  1 1/12/94  and  at  sites  C2  and  C4  at  2300  hr,  1 1/10/94. 


During  peak  current,  due  to  significant  vertical  turbulent  diffusion  and  strong  upward 
mixing  over  the  lutocline  (caused  by  internal  wave  breaking  and  high  rate  of  erosion  of 
freshly  deposited  sediment  during  the  previous  slack),  SSC  over  the  water  column  usually 
increases,  the  elevation  of  lutocline  layer  rises,  and  both  fluid  mud  and  lutocline  layers 
become  comparatively  thicker  (Odd  and  Rodger,  1986;  Wolanski,  et  al.,  1988).  In  the 
Jiaojiang,  as  noted,  although  the  ebb  peak  is  weaker  than  flood,  it  lasts  longer  than  flood.  As 
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a resiilt,  during  peak  ebb,  diffusion,  vertieal  mixing  over  the  lutocline  and  bottom  erosion  are 
also  comparatively  high.  Consequently,  both  during  peak  ebb  and  peak  flood,  SSC  exhibits 
quantitatively  the  same  behavior  as  noted,  except  near  the  bottom  during  peak  flood  due  to 
the  cumulative  effects  of  dilution,  stronger  entrainment  over  the  lutocline  and  mixing  within 
the  fluid  mud  layer. 

6.4.4  Vertical  Profiles  of  Velocity 

Typical  vertical  profiles  of  velocity  over  a tidal  cycle  are  shown  in  Figures  6.1 1-6.14. 
Many  of  the  simulated  velocity  profiles  are  in  reasonable  agreement  with  the  observations, 
although  the  simulations  for  spring  tide  show  better  agreement  than  at  neap.  It  is  also 
observed  that  at  peak  flood,  the  velocity  over  the  entire  water  column  was  greater  than  that 
at  peak  ebb.  Flood  peak  near  the  bottom  ranges  0.8-1. 0 m s at  spring  tide  and  0.5-0.8 

m s at  neap.  Ebb  peak  near  the  bottom  ranges  0.6-0.8  m s at  spring  tide  and  0.4-0.65 
m s'*  at  neap. 

Noteworthy  differences  between  simulations  and  observations  are  as  follows;  during 
the  flood  at  neap  tide,  the  model  over-predicted  some  profiles  near  the  bottom,  e.g.,  at  2 hr 
and  3 hr  in  Figures  6.12  and  6.14,  respectively.  During  ebb,  the  model  over-predicted  some 
profiles  near  the  surface  at  C2,  e.g.,  at  8 hr  and  9 hr  in  Figure  6.12  and  under-predicted  some 
profiles  near  the  surface  at  C4,  e.g.,  at  9 hr  and  10  hr  in  Figure  6.14.  These  differences 
between  simulations  and  observations  are  perhaps  caused  by  following  factors: 

(i)  Use  of  Approximate  Stratification  Function:  The  stratification  function  [Eq. 
(3.10)]  used  in  the  model  can  be  expected  to  affect  the  bottom  drag  coefficient,  [Eq.  (3.9) 
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Figure  6.11.  Velocity  profiles  at  site  C2  during  a spring  tide.  Figure  6.12.  Velocity  profiles  at  site  C2  during  a neap  tide. 

Solid  lines  are  simulations  and  open  circles  represent  data  Solid  lines  are  simulations  and  open  circles  represent  data 

obtained  during  1900  hr,  1 1/05/94  to  0600  hr,  1 1/06/94.  obtained  during  1 100  hr  to  2100  hr,  1 1/10/94.  Positive 

Positive  values  signify  flood  and  negative  denote  ebb.  values  signify  flood  and  negative  denote  ebb. 
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Solid  lines  are  simulations  and  open  circles  represent  data  Solid  lines  are  simulations  and  open  circles  represent  data 

obtained  during  1900  hr,  1 1/05/94  to  0500  hr,  1 1/06/94.  obtained  during  1 100  hr  to  2100  hr,  1 1/10/94.  Positive 

Positive  values  signify  flood  and  negative  denote  ebb.  values  signify  flood  and  negative  denote  ebb. 
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and  Table  6.1],  and  consequently  the  near-bottom  velocity  near  the  bottom.  As  a result,  an 
over-prediction  of  the  stratification  effect  can  imply  reduced  and  enhanced  velocity  near 

the  bottom. 

(2)  Coriolis  Effect:  The  model  omitted  the  Coriolis  effect  due  to  the  high  Rossby 
number  within  the  modeled  domain,  R^-U/f  5=12.0»1,  where  the  Coriolis  parameter 

/=2Qsin(t),=  7.29x10'^  s the  angular  velocity  of  the  earth’s  rotation  Q =7.29x10"^  s’, 

the  latitude  (j), -30° , the  characteristic  mean  velocity  C/=1.0  m s and  the  mean  width  of 

the  estuary  b =1.2  km  (Section  4.1).  Note  that  signifies  the  effect  of  flow  inertia  relative 

to  Coriolis  acceleration.  Thus,  a high  R^  implies  comparatively  low  Coriolis  effect.  Note 

however,  that  Taizhou  Bay  outside  the  estuary  (Figure  4.1)  is  a much  wider  water  body  with 
a mean  width  of  about  6 km  and  R^  of  about  2.0.  Thus,  one  may  expect  that  Coriolis 

acceleration  would  have  a greater  effect  on  the  tidal  current  in  the  bay,  which  in  turn  perhaps 
has  an  effect  on  the  tidal  currents  within  the  estuary.  However,  in  that  context  the  following 
analysis  should  be  noted. 

To  examine  the  Coriolis  effect  fiirther,  the  work  of  Proudman  (1925)  is  introduced 
here.  Proudman  calculated  the  vertically-averaged  transverse  and  longitudinal  tidal  currents, 
V and  U,  across  a channel  under  the  Coriolis  effect  for  the  special  case  of  the  angular 

frequency  of  the  n'^  tidal  constituent,  which  corresponds  to  the  diurnal  tide,  S2,  at  latitude 


30°.  A parabolic  section  is  assumed  for  the  channel  as 
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(6.1) 

where  h is  the  mean  water  depth  below  MSL,  is  the  maximum  water  depth  at  the  channel 

center  and>^  is  the  transverse  Cartesian  coordinate  originating  from  the  channel  center.  The 
relevant  formulas  are  as  follows: 

Non-dimensional  transverse  current: 


O)  pkB 

JLhV=  ^ 


SK  Ik^'^ 


2(0^  . , - 
sK 


, (6.2) 


sinhA(5  +y)  +[A:  -y  ^) -k(B  +;/)]coshA(5  +y) 


Non-dimensional  longitudinal  current: 


-^hU=-  ^ 


SK  2k^'^ 


2o)^  , , - 

l^—^+k\W-y^)-k(B-y) 

sK 


(6.3) 


sinhA(5  +>^)  -[A:  -y  +k{B  +>;)]coshA(5  +y) 


with  the  tidal  wave  number  k=(jiy\[^,  where  h is  the  mean  water  depth  of  the  parabolic 


flow  cross-section,  i.e.,  h=2h^3 . Figure  6.15  shows  the  distributions  of  the  transverse  and 
longitudinal  currents  across  channel  in  the  Jiaojiang  and  also  Taizhou  Bay.  h^=5  m and 
(0^=7.27x10'^  rad  s'*  (diurnal  tide)  were  taken  both  for  the  estuary  and  the  bay,  and  b^= 1.2 


km  and  6 km  were  taken  in  the  estuary  and  the  bay,  respectively.  It  is  observed  that  even  in 
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the  Taizhou  Bay  the  Coriolis  effect  is  almost  negligible.  The  circulation  in  the  bay  due  to  this 
effect  will  rotate  in  the  counterclockwise  sense.  This  will  lead  to  a clockwise  circulation 
within  the  estuary  due  to  the  vorticity  transfer  by  shear,  consistent  with  the  sense  of  rotation 
observed  (Section  6.4.2). 


y/B 


(a) 


(b) 


Figure  6.15.  Distributions  of  transverse  (a)  and  longitudinal  (b)  currents  across  the 
channel  due  to  the  Coriolis  effect,  viewed  in  the  direction  of  tidal  wave  propagation. 
Solid  lines  are  currents  in  Taizhou  Bay  and  dashed  lines  in  the  Jiaojiang. 


(3)  Gravitational  Circulation  due  to  Salinity:  Previous  observations  in  the  Jiaojiang 
have  showed  that  salinity  is  typically  well-mixed  vertically  during  the  entire  tidal  cycle 
(Zhou,  1986;  Li,  et  al.,  1993;  Dong,  et  al.,  1997).  Hence,  salinity-induced  transport  was  not 
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considered  in  the  present  modeling  effort.  However,  based  on  their  simulations  of  salinity 
distribution  in  the  Jiaojiang  with  and  without  sediment-induced  buoyancy  effects,  Guan,  et 
al.  (1998)  observed  that  suspended  sediment  could  measurably  modify  the  vertical 
distribution  of  salinity  due  to  the  fact  that  suspended  sediment  causes  the  currents  to  be 
strongly  sheared  along  the  vertical  direction.  They  predicted  that  under  sediment-induced 
buoyancy  effects,  salinity  stratification  locally  increased  during  flood  tide,  thus  in  turn 
modifying  the  vertical  profile  of  velocity. 

Gravitational  circulation  due  to  salinity  stratification  can  be  approximately  expressed 
in  a formulation  including  river  inflow  and  surface  wind  stress  as  follows  (Hansen  and 
Rattray,  1965): 

^ ^ nn  R.CI 

^=_(l  -o2)+l(l  +40+30^)+^— (1  -9o2-8o3)  (6.4) 

where  the  first  term  on  the  right  hand  side  is  due  to  the  effect  of  river  inflow,  the  second  term 
is  due  to  the  effect  of  wind,  the  third  term  is  due  to  salinity-induced  density  gradient,  is 

the  gravitational  circulation  and  Uj.  is  the  vertical-mean  inflow  velocity.  Further,  T is  the 

dimensionless  wind  stress,  equal  to  is  the  river  inflow  rate,  Ra  is  the 

estuarine  Rayleigh  number,  defined  as  the  ratio  of  free  convection  to  diffusion  (Hansen  and 
Rattray,  1965),  i.e.,  equal  to  is  a proportionality  coefficient  between 

fluid  density  and  salinity  (=0.78x  10’^),  is  the  salinity  at  the  mouth  of  the  estuary,  and 

is  a constant  in  the  relationship  of  salinity  distribution  introduced  as  follows 
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X+Z(o)  .. 

Sq  ^ (6-5) 

Here  X is  the  dimensionless  horizontal  coordinate,  i.e.,  X=R^/BHA^,  and  Z is  the  vertical 
distribution  flmction  of  salinity. 

Figure  6.16  shows  profiles  of  for  different  a^Ra  without  considering  the  effect 

of  surface  wind  stress.  This  type  of  gravitational  circulation  can  cause  ebb  current  to 
increase  near  the  surface  and  decrease  near  the  bottom,  and  vice  the  versa  during  flood.  The 
observations  at  C4  (Figure  6.13  and  14)  were  characteristically  consistent  with  this  trend.  It 
should  be  noted  that  in  light  of  the  likely  interaction  between  salinity  and  SSC  in  the 
Jiaojiang,  Hansen  and  Rattray’s  analysis  based  on  salinity  effect  alone  must  be  interpreted 
more  broadly  for  this  estuary  with  regard  to  the  likely  combined  effects  of  the  salinity  and 
SSC  in  inducing  the  gravitational  circulation  shown  in  Figure  6.16 

(4)  Geomorphologic  Effects:  Due  to  non-linear  effects  of  bottom  friction  and  of 
advection,  tidal  residual  currents  usually  exist  in  the  vicinity  of  topographies  such  as  basins, 
headlands,  sand  ridges,  etc.  (Zimmerman,  1981).  Noteworthy  features  in  the  Jiaojiang  are  the 
headlands  located  at  the  mouth  of  Taizhou  Bay  (Figure  4.1).  It  is  likely  that  tidal  residual 
currents  induced  by  these  headlands  have  effects  on  the  residual  currents  within  the  estuary 
itself  As  shown  in  Figure  6.17  (Zimmerman,  1981),  along  a coastal  promontory  tidal  flow 
tends  to  accelerate  as  well  as  decelerate,  with  a maximum  near  the  headland.  Concurrently, 
a frictional  boundary  layer  develops  with  diminishing  current  velocity  towards  the  coast. 
Thus,  vorticity  is  generated  along  the  coast  with  orientation  alternating  between  the  flood  and 
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ebb.  As  the  vorticity  is  largest  near  the  headland,  during  flood  there  is  a net  flux  of 
counterclockwise  vorticity  into  the  left  closed  curve  and  out  of  the  right  one.  During  the  ebb 
the  situation  does  not  reverse  because  of  a net  flux  of  clockwise  vorticity  out  of  the  left 
quadrilateral  and  into  the  right  quadrilateral,  which  is  equivalent  to  the  situation  during  the 
flood.  Hence  there  is  a net  flux  of  counterclockwise  vorticity  into  the  left  quadrilateral  and 
of  clockwise  vorticity  into  the  right  quadrilateral,  giving  rise  to  a vortex  pair  of  opposite 
signs  at  either  sides  of  the  headland.  In  other  words,  residual  circulation  exists  at  both  sides 
of  the  headland. 


-15  -10  -5  0 5 10  15  20 


u,/U, 

Figure  6.16.  Profiles  of  gravitational  circulation  without 
wind  stress,  after  Hansen  and  Rattray  (1965). 
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Thus,  a clockwise  circulation  can  be  expected  in  the  north  region  of  Taizhou  Bay  and 
a counterclockwise  one  in  its  southern  region.  These  circulations  possibly  have  extended 
effects  on  residual  currents  within  the  Jiaojiang.  Unfortunately,  Taizhou  Bay  is  beyond  the 
modeled  domain,  and  there  are  no  observations  available  to  back-up  this  assertion. 


Figure  6.17.  Tidal  currents,  vorticities  and  residual  circulation  in  the  neighborhood  of  a 
headland.  Currents  are  signified  by  solid  arrows  for  flood  and  dashed  arrows  for  ebb. 
Currents  are  largest  near  the  headland  and  decrease  towards  the  shoreline.  Vorticity 
therefore  has  a maximum  near  the  headland.  Vorticities  generated  by  side-wall  fnction 
are  shown  by  solid  circles  for  flood  and  dashed  circles  for  ebb.  Vorticities  have  highest 
strength  near  the  headland  and  diminish  away  from  it  (after  Zimmerman,  1981). 


6.4.5  Vertical  Profiles  of  SSC 

Figures  6.18-6.21  show  vertical  profiles  of  SSC  simulated  and  observed  at  sites  C2 
and  C4  over  one  tidal  cycle  during  spring  and  neap  tides,  respectively.  The  following  features 
of  these  SSC  profiles  are  noteworthy: 

(1)  SSC  changes  in  the  following  way:  near  the  bottom  during  ebb,  SSC  continues 
to  increase  until  flow  deceleration  occurs,  then  decreases.  During  flood  it  either  continues 
to  decrease,  or  increases  somewhat  during  the  accelerating  phase  of  flow,  then  decreases. 
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lines  are  simulations  and  open  circles  represent  field  data  lines  are  simulations  and  open  circles  represent  field  data 

obtained  during  1900  hr,  1 1/05/94  to  0600  hr,  1 1/06/94.  obtained  during  1 100  hr  to  2100  hr,  1 1/10/94. 
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Near  the  surface,  it  increases  during  accelerating  flow  and  decreases  during  the  decelerating 
flow  both  for  flood  and  ebb. 

(2)  During  peak  currents,  the  lutocline  layer  (defined  in  Section  2.1)  rises  and 
becomes  thicker,  with  gentler  upper  and  lower  gradients.  During  slack  water  it  fall  and 
becomes  thinner,  with  steeper  upper  and  lower  gradients. 

(3)  The  lutocline  layer  is  lower  and  thinner  during  neap  tide  than  during  spring. 

(4)  Fluid  mud  and  lutocline  layers  occur  over  the  entire  tidal  cycle  except  during  the 
period  around  high  water  slack. 

This  behavior  of  SSC  reflects  the  cumulative  effects  of  floe  settling,  tidal  asymmetry, 
lutocline  entrainment,  turbulence  damping  in  fluid  mud  layer  and  the  large  ratio  of  tidal 
range  to  water  depth  (see  also  Section  6.4.3). 

Differences  between  the  simulations  and  observations  of  SSC  are  as  follows:  during 
ebb  tide,  the  model  over-predicts  SSC  for  some  profiles  near  the  bottom,  e.g.,  5 hr,  7 hr  and 
1 1 hr  in  Figure  6. 1 9,  and  under-predicts  near  the  surface,  e.g.,  1 0 hr  and  1 1 hr  in  Figure  6.18 
and  8 hr,  9 hr  and  10  hr  in  Figure  6.20.  During  flood  tide,  the  model  under-predicts  SSC  in 
some  profiles  over  the  entire  water  column,  e.g.,  2 hr  and  3 hr  in  Figures  6.18  and  6.20, 
respectively.  These  differences  likely  arise  due  to  following  modeling  approximations  as  well 
as  observation  errors: 

(1)  Approximate  Open  Boundary  Conditions  for  SSC:  SSC  at  the  open  boundaries 
were  measured  30  hr  before  and  after  those  at  sites  C2  and  C4  during  spring  and  neap  tides, 
respectively,  and  were  taken  as  uniform  over  the  cross-sections  of  the  open  boundaries  (see 
also  in  Section  6.3.3). 
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In  order  to  demonstrate  the  lateral  non-uniformity  of  SSC  in  the  Jiaojiang,  model 
simulations  and  observations  of  vertically-averaged  SSC  at  the  flow  section  including  sites 
C2  and  C4  are  plotted  in  Figure  6.22.  From  this  figure  together  with  the  tidal  variation  of 
velocity  [Figures  6.6.6(a)-6.9(a)],  it  is  found  that  during  flood,  SSC  at  C4  was  greater  than 
at  C2,  and  during  ebb  it  gradually  increased  at  C2  and  ultimately  became  greater  than  at  C4. 
A difference  (in  the  vertically-averaged  SSC)  of  about  0-7  kg  m occurred  between  C2  and 

C4,  which  accounts  for  about  0%-70%  of  the  vertically-averaged  values  of  SSC  at  these  sites. 
Because  uniform  SSC  over  the  cross-sections  of  the  open  boundaries  were  employed  in 
modeling,  the  model  generated  a more  uniform  SSC  across  the  estuary  within  the  modeled 
domain  than  observations  (Figure  6.22).  Table  6.3  includes  typical  values  from  Figure  6.22 
at  4 hr,  where  two  exhibit  significant  SSC  non-uniformity  and  the  other  two  do  not.  It  is 
evident  that  better  comparisons  of  vertically-averaged  SSC  between  simulations  and 
observations  resulted  at  times  of  comparatively  minor  non-uniformity  than  when  significant 
lateral  non-uniformity  occurred. 


Table  6.3.  Vertically-averaged  SSC  during  minor  and  significant 
non-uniformity  of  SSC  across  the  flow  cross-section  (Unit:  kg  m '^) 


Method 

Site 

Minor  non-uniformity 

Significant  non-imiformity 

Spring  tide 

Neap  tide 

Sprin 

g tide 

Neap  tide 

5hr 

lOhr 

5hr 

8hr 

2hr 

15hr 

2hr 

lOhr 

Observation 

C2 

5.50 

12.50 

4.00 

6.00 

11.50 

6.75 

6.75 

13.20 

C4 

6.50 

12.90 

6.00 

7.00 

16.25 

14.00 

12.10 

9.00 

Simulation 

C2 

5.70 

13.25 

4.00 

7.00 

11.50 

10.75 

13.50 

12.50 

C4 

5.70 

15.00 

3.60 

7.50 

9.50 

11.75 

14.00 

15.50 

160 


Figure  6.22.  Time  series  of  vertically-averaged  SSC.  Dark  circles  are  from 
site  C2,  open  circles  represent  site  C4,  solid  lines  signify  simulations  at  C2, 
dashed  lines  are  simulations  at  the  center  of  the  flow  section  containing  C2 
and  C4  and  dotted  lines  denote  simulations  at  C4. 


In  order  to  demonstrate  the  significance  of  SSC  non-uniformity,  a model  test  using 
non-uniform  boundary  conditions  for  SSC  was  carried  out  for  a spring  tide.  The  simulated 
results  are  shown  in  Figure  6.23.  Also  shown  in  this  figure  are  results  using  uniform  SSC 
boundary  conditions  (Figure  6.22).  Non-uniform  SSC  at  the  boundaries  was  generated  using 
a linear  relationship  as  follows:  (1)  the  mean  value  of  SSC  over  each  open  boundary  was 
taken  as  the  observation  at  site  Cl  or  C3;  (2)  the  slope  of  the  lateral  variation  of  SSC  was 
taken  as  the  ratio  of  vertically-averaged  SSC  observed  between  sites  C2  and  C4  (Figure 
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6.22).  From  Figure  6.23,  it  is  observed  that  considerable  differences  in  the  simulated  SSC 
were  generated  during  peak  ebb.  This  suggests,  at  least  qualitatively,  that  non-uniformity  of 
SSC  at  the  open  boundaries  may  have  had  a noteworthy  effect  on  simulations  within  the 
estuary. 


Time  starting  2000  hr,  1 1/10/94. 

Figure  6.23.  Comparision  of  simulated  vertically-averaged  SSC  using  uniform 
and  nonuniform  boundary  conditions  of  SSC  during  a spring  tide.  Solid  lines 
signify  simulations  at  site  C2,  dashed  lines  are  that  at  site  C4  and  dark  and 
open  circles  represent  that  using  uniform  boundary  conditions  of  SSC. 


(2)  Approximate  Parameters:  The  selection  of  modeling  parameters  given  in  Table 
6. 1 was  a difficult  task.  Although  most  parameters  were  determined  optimally  from  previous 
works,  data  analysis  and  model  calibration,  parametric  selection  can  in  general  cause 
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measurable  differences  between  simulations  and  observations  due  to  the  approximate  nature 
of  the  parameters.  Noted  in  the  following  are  those  parameters  that  are  sediment-dependent 
and  highly  influential  in  controlling  the  sensitivity  of  the  simulations  of  SSC; 

Parameters  Determined  bv  Model  Calibrations:  Model  calibrations  were  carried  out 
based  on  comparisons  of  SSC  between  simulations  and  observations.  Because  all  parameters 
affected  each  other  on  an  inter-dependent  way  in  the  calibration  tests,  it  was  at  times  difficult 
to  identify  the  desired  value  of  a certain  parameter.  In  this  regard,  the  most  noteworthy  ones 
are;  the  efficiency  coefficient  of  turbulence  damping,  d^,  in  the  vertical  mass  diffusion 

equation  (3 .20)  (see  also  in  Table  6. 1 ),  and  in  the  consolidation 

rate  formula  (2.23)  (see  also  in  Table  6.1)  and  C,  Pj  and  in  the  vertical  distribution 

formula  (2.24)  of  dry  sediment  concentration  for  a fully-consolidated  bed  (Table  6.1).  There 
were  no  direct  experimental  data  available  for  determining  these  parameters  for  the  Jiaojiang. 

Bottom  Erosion  Rate  Constant:  The  bottom  erosion  rate  constant,  was  obtained 

from  analysis  of  field  data  (Section  4.3.6  and  Table  6.1).  Data  analysis  was  carried  out  based 
on  three  basic  assumptions  (Section  4.3.6).  It  is  known  that  the  second  and  third  assumptions 
seem  reasonable  but  the  first  my  not  be  so.  As  stated  in  Section  6.4.3,  advection  is  a major 
factor  affecting  sediment  transport  in  the  Jiaojiang,  and  consequently  has  a significant  effect 
on  SSC.  Thus  it  can  not  be  ignored.  In  fact,  the  scatter  of  data  points  in  Figure  4.29 
demonstrates  the  effect  of  advection.  Thus,  the  correct  erosion  rate  constant  can  only  be 
determined  through  erosion  experiments  in  laboratory  using  the  same  mud  as  in  the 
Jiaojiang. 
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(3)  Observation  Errors:  As  noted  in  Chapter  4,  a ship-borne  turbidimeter  was  used 
for  measuring  the  SSC.  The  ship  was  anchored  by  a mooring  chain  with  a length  = 10//=50 
m.  Thus,  the  ship  could  turn  along  with  tidal  currents  and  change  its  orientation  within  the 
circle  defined  by  the  chain.  Given  the  laterally  non-uniform  distribution  of  SSC  in  the 
Jiaojiang  as  described  above  and  in  Section  6.4.3,  changing  ship  position  could  have  led  to 
a degree  of  error  in  the  observed  SSC. 

6.4.6  Lutocline  Laver 

The  lutocline  layer  was  identified  by  its  definition  in  Section  2. 1 , i.e.,  the  near-bottom 
layer  between  the  upper  elevation  and  lower  elevation  of  the  maximum  vertical  gradient  of 
SSC  (Figure  2.1).  Figure  6.24  shows  typical  vertical  distributions  of  the  vertical  gradient  of 
SSC.  The  lutocline  layer  is  identified  by  the  zone  over  which  the  SSC  gradient  is 
comparatively  high. 

Time  variations  of  this  layer,  both  simulated  and  observed  at  sites  C2  and  C4  during 
spring  and  neap  tides,  are  shown  in  Figures  6.24-6.27.  It  is  seen  that  both  the  measured 
pattern  and  elevation  of  the  lutocline  layer  are  reproduced  approximately  by  the  model. 
Except  for  periods  of  about  1-2  hr  around  high  water  slack,  a lutocline  layer  with  a thickness 
of  about  1-3  m (Figure  6.29)  and  an  upper  elevation  of  1-4  m (Figure  6.30)  consistently 
occurred,  irrespective  of  whether  the  tide  was  spring  or  neap. 

From  Figures  6.25-6.28  and  Table  6.4,  combined  with  Figures  6.6(a)-6.9(a),  the 
following  characteristics  of  the  lutocline  layer  can  be  gleaned: 

(1)  The  lutocline  layer  was  at  higher  elevation  during  peak  flows  than  during  slack 
water,  with  the  highest  elevations  during  peak  floods. 
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Figure  6.24.  Vertical  gradient  of  SSC  as  a function  of  elevation  during 
a neap  tide.  Solid  lines  are  simulations  and  open  circles  represent 
field  data  obtained  during  1 100  hr  to  2100  hr,  1 1/10/94. 


(2)  The  lutocline  layer  was  thinner  during  slack  water  than  during  peak  flow,  and  also 
thinner  during  neap  tide  than  during  spring. 

(3)  The  lutocline  layer  rose  during  period  of  flow  acceleration  and  fell  during  flow 
decelerating  periods. 

From  comparisons  of  the  thickness  of  the  lutocline  layer  between  simulations  and 
measurements  (Figure  6.29),  it  is  observed  that  the  predicted  thickness  of  the  lutocline  layer 
is  in  better  agreement  with  measurement  at  neap  tide  than  at  spring  (Figure  6.29).  Also,  there 
is  a degree  of  over-prediction  of  its  upper  limit  at  neap  tide  and  under-prediction  at  spring 
(Figure  6.30). 


Elevation  above  bed  (m)  Elevation  above  bed  (m) 
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Figure  6.24.  Simulated  (solid  lines)  and  measured  (dashed  lines) 
tidal  variation  of  the  lutocline  layer  at  site  C2  during  a spring 
tide  from  1800  hr,  1 1/05/94  to  1700  hr,  1 1/06/94. 


Figure  6.25.  Simulated  (solid  lines)  and  measured  (dashed  lines) 
tidal  variation  of  the  lutocline  layer  at  site  C2  during  a neap 
tide  from  1000  hr,  1 1/10/94  to  0900  hr,  1 1/1 1/94. 


Elevation  above  bed  (m)  Elevation  above  bed  (m) 
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Figure  6.26.  Simulated  (solid  lines)  and  measured  (dashed  lines) 
tidal  variation  of  the  lutocline  layer  at  site  C4  during  a spring 
tide  from  1800  hr,  1 1/05/94  to  1700  hr,  1 1/06/94. 


Figure  6.27.  Simulated  (solid  lines)  and  measured  (dashed  lines) 
tidal  variation  of  the  lutocline  layer  at  site  C4  during  a neap 
tide  from  1000  hr,  1 1/10/94  to  0900  hr,  1 1/1 1/94. 
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Figure  6.28.  Lutocline  layer  thickness:  simulation  (6,J  and  measurement  (6,J. 
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Table  6.4.  Average  thickness  and  upper  elevation 
__ofTutocline  layer  at  different  times  (Unit:  m) 


Item 

Tidal  range 

Observation 

Simulation 

Spring  tide 

Nea 

3 tide 

Spring  tide 

Neap  tide 

C2 

C4 

C2 

C4 

C2 

C4 

C2 

C4 

Lutocline 

layer 

thickness 

Peak  flood 

1.83 

1.61 

1.17 

1.02 

1.98 

1.83 

1.46 

1.42 

Peak  ebb 

1.46 

1.54 

1.10 

1.02 

1.61 

1.24 

1.10 

1.17 

High  slack 

0.00 

0.00 

0.58 

0.88 

0.00 

0.00 

0.44 

1.02 

Low  slack 

1.32 

1.17 

1.02 

0.80 

1.17 

1.02 

1.02 

0.88 

Tidal  mean 

1.27 

1.30 

0.97 

0.97 

1.29 

1.43 

0.99 

1.33 

Upper 
elevation 
of  lutocline 
layer 

Peak  flood 

3.07 

3.80 

2.63 

3.07 

3.29 

3.59 

3.15 

3.50 

Peak  ebb 

3.15 

3.80 

1.98 

1.54 

3.15 

3.07 

2.34 

2.56 

High  slack 

0.00 

0.00 

1.10 

1.17 

0.00 

0.00 

0.59 

1.17 

Low  slack 

1.90 

2.20 

1.61 

1.61 

2.05 

1.76 

1.61 

2.05 

Tidal  mean 

2.41 

3.31 

1.64 

1.80 

2.44 

2.41 

2.09 

2.51 

6.4.7  Flow-SSC  Hysteresis 

Due  to  time-lag  effects  associated  with  settling,  diffusion,  bed  erosion,  entrainment 
and  consolidation,  SSC  in  estuaries  is  usually  higher  during  decreasing  currents  than  when 
currents  are  increasing  (Postma,  1967;  Dyer  and  Evans,  1989).  As  a result,  when  the  SSC  at 
a certain  elevation  is  plotted  against  the  bottom  shear  stress,  this  hysteresis  becomes  visually 
evident  (Costa  and  Mehta,  1990). 

Simulated  and  measured  hysteresis  loops  of  SSC  for  sites  C2  and  C4  are  shown  in 
Figures  6.31-6.38.  These  are  at  1 m above  the  bottom  and  1 m below  the  instantaneous 
surface  during  spring  and  neap  tides. 
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Figure  6.30.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at  1 m 
above  bottom  at  C2  during  a spring  tide  from  1800  hr,  1 1/05/94  to  0500  hr,  1 1/06/94. 


(Pa) 


Figure  6.31.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at  1 m 
below  surface  at  C2  during  a spring  tide  from  1800  hr,  1 1/05/94  to  0500  hr,  1 1/06/94. 
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Figure  6.32.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at 
1 m above  bottom  at  C2  during  a neap  tide  from  1000  hr  to  2100  hr,  1 1/10/94. 


(Pa) 

Figure  6.33.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at 
1 m below  surface  at  C2  during  a neap  tide  from  1000  hr  to  2100  hr,  1 1/10/94. 
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Figure  6.34.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at  1 m 
above  bottom  at  C4  during  a spring  tide  from  1800  hr,  1 1/05/94  to  0500  hr,  1 1/06/94. 


•^4  (Pa) 


Figure  6.35.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at  1 m 
below  surface  at  C4  during  a spring  tide  from  1800  hr,  1 1/05/94  to  0500  hr,  1 1/06/94. 
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Figure  6.36.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at 
1 m above  bottom  at  C4  during  a neap  tide  from  1000  hr  to  2100  hr,  1 1/10/94. 


X*  (Pa) 

Figure  6.37.  Simulated  (solid  line)  and  measured  (dashed  line)  hysteresis  loops  at 
1 m below  surface  at  C4  during  a neap  tide  from  1000  hr  to  2100  hr,  1 1/10/94. 
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During  ebb  (Figures  6.31-6.38),  the  simulated  and  measured  loops  near  the  bottom 
and  near  the  surface  follow  the  same  trend.  During  accelerating  flow  SSC  increases  gradually 
(by  resuspension).  As  the  flow  begins  to  decelerate,  at  first  SSC  increases  dramatically  until 
it  reaches  a maximum,  then  starts  to  decrease  (due  to  deposition).  This  behavior  results  in 
the  loops  rotation  in  the  clockwise  sense. 

During  flood,  the  loops  near  the  bottom  and  near  the  surface  follow  different  patterns. 
Near  the  bottom  (Figures  6.31,  6.33,  6.35  and  6.37),  SSC  either  decreases  over  all  stages,  or 
increases  somewhat  during  accelerating  flow,  then  decreases.  This  causes  the  loops  to  rotate 
clockwise.  Near  the  surface  (Figures  6.32,  6.34,  6.36  and  6.38),  SSC  increases  during 
accelerating  flows  due  to  vertical  diffusion  and  mixing  and  reaches  a maximum  at  the  end 
of  the  period  of  acceleration.  Subsequently,  during  decelerating  flows,  SSC  begins  to 
decrease  due  to  settling  and  dilution.  This  causes  the  loops  rotating  in  the  counterclockwise 
sense. 

Thus,  during  flood,  loop-reversals  appeared  near  the  bottom  in  the  Jiaojiang.  Costa 
and  Mehta  (1990)  observed  such  loop-reversals  in  Hangzhou  Bay,  China,  during  the 
transition  from  accelerating  to  decelerating  flow.  They  showed  that  reversal  was  induced  by 
the  presence  of  the  lutocline  at  that  level.  During  accelerating  flow  the  lutocline  forms,  and 
entrainment  over  it  occurs  concurrently.  Due  to  the  entrainment  lag,  SSC  near  the  bottom 
increases  and  reaches  a maximum  during  accelerating  flow,  then  decreases.  This  leads  to  a 
loop-reversal. 

In  the  Jiaojiang,  more  significant  loop-reversals  occurred  near  the  bottom  during  the 
entire  flood  tide  period  (Figures  6.31,  6.33,  6.35  and  6.37)  with  the  SSC  either  decreasing 
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over  this  period,  or  increasing  somewhat  during  aecelerating  flow,  then  decreasing.  As  noted 
in  section  6.4.3,  this  behavior  is  believed  to  be  caused  by  the  effects  of  dilution,  vertical 
diffusion  as  well  as  the  presence  of  the  lutocline. 

6-4.8  Effect  of  Turbulence  Damping  on  SSC  and  Lutocline  Formation 

With  reference  to  Section  6.4.3,  in  order  to  further  examine  the  effect  of  turbulence 
damping  on  SSC  and  the  lutocline  formation,  the  following  numerical  tests  were  carried  out: 

(1)  Test  1 : Considering  turbulence  damping  and  taking  d^=Q.lS ; 

(2)  Test  2:  Neglecting  turbulence  damping,  i.e.,  taking 

(3)  Tgst  3:  Neglecting  the  effects  of  cohesion  and  interactions  between  floes,  i.e., 
taking  Ri^^Ri^^O  and  d^^Q.lS. 

In  the  above  tests,  all  other  parameters  were  taken  to  be  the  same  as  in  previous 
simulations.  Figure  6.39  shows  the  results  of  these  tests  along  with  the  observations  at  C4 
during  a neap  tide.  It  is  seen  that  without  considering  turbulence  damping  (Test  2),  the  model 
resulted  in  relatively  uniform  profiles  of  SSC  and  consequently  less  distinct  lutoclines  (2-5 
hrs  and  9-1 1 hrs  in  Figure  6.39).  If  compared  with  the  observations,  it  is  seen  that  Test  2 had 
the  worst  results.  Test  1 had  the  best  with  Test  3 in-between.  However,  there  were  smaller 
differences  among  these  tests  during  low  slack  water  (1  hr  in  Figure  6.39)  and  high  slack 
water  (7  and  8 hrs  in  Figure  6.39).  This  implies  that  during  high  and  low  slack  waters,  the 
lutocline  is  mainly  governed  by  floe  settling.  It  is  also  seen  that  considerable  effect  of 
turbulence  damping  appeared  near  the  bottom  (2-5  hrs  and  9-1 1 hrs  in  Figure  6.39),  because 
of  damping  induced  by  the  relatively  high  SSC  there.  Near  the  surfaee,  the  effect  of  damping 
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during  flood  (2-5  hrs  in  Figure  6.39)  is  greater  than  during  ebb  (9-1 1 hrs  in  Figure  6.39),  due 
to  the  relatively  higher  SSC  near  the  surface  during  flood.  Without  considering  cohesion  and 
interactions  between  floes  (Test  3),  the  simulated  results  are  very  close  to  Test  1 and  the 
differences  are  mainly  near  the  bottom.  This  supports  the  conclusion  in  Chapter  5 that 
cohesion  and  interactions  between  floes  govern  turbulence  damping  only  at  high  SSC  (>40- 
50  kg  m ^),  in  the  hindered  settling  range. 
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Figure  6.39.  Modeling  SSC  profiles  at  site  C4  during  a neap  tide.  Open 
circles  represent  field  data  from  1 100  hr  to  2100  hr,  1 1/10/94,  solid  lines 
are  simulations  with  d2=0.75,  dashed  lines  signify  simulations  with  d2  = 

0.75  and  Rig=Ri^=0,  and  dotted  lines  represent  simulations  with  d2=0. 


CHAPTER  7 

SUMMARY  AND  CONCLUSIONS 
7.1  Summary 

Lutoclines  are  common  features  in  estuaries  vvith  high-load  fine  sediments.  They  are 
significantly  associated  with  interfacial  waves  and  turbulenee  damping  in  fluid  mud,  and  tend 
to  diminish  vertical  mixing  under  tidal  forcing.  In  order  to  better  understand  these  dynamical 
features,  the  following  studies  were  carried  out: 

1 . Turbulenee  Damping:  Following  Rossby  and  Montgomery  (1935),  assuming  the 
sum  of  the  kinetic,  potential  and  dissipated  energy  in  a stratified,  cohesive  flow  to  be  equal 
to  that  in  a homogeneous,  non-cohesive  flow  at  identical  shear  rates,  turbulence  damping  in 
fluid  mud  was  phenomenologically  examined.  Accordingly,  sediment  settling  flux,  cohesion, 
interaetions  between  floes  and  sediment-induced  stratifieation  were  quantified  by  the 
respective  Richardson  numbers  Ri^,  Ri^,  Ri^  and  Ri.  The  resulting  formulations  for  the 

turbulent  mixing  length  were  examined  using  observations  in  the  Jiaojiang  estuary.  The 
corresponding  expressions  for  the  vertical  momentum  and  mass  diffusion  coefficients  were 
incorporated  in  the  developed  numerical  model  codes  for  flow  and  sediment  transport. 

2.  Internal  Waves:  ASSM  data  from  the  Jiaojiang  estuary  were  used  to  examine  the 
height,  angular  frequency,  celerity  and  length  of  internal  waves  at  the  lutocline.  Their 
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variation  with  the  global  Richardson  number  was  examined  in  the  light  of  results  of  previous 
studies  on  the  lutocline  and  other  pycnoclines. 

3-  Lutocline  Response  to  Tidal  Forcing:  Three-dimensional,  finite-difference 
hydrodynamic  and  sediment  transport  codes,  COHYD-UF  and  COSED-UF,  were  developed, 
incorporating  the  latest  unit  process  models  for  fine  sediment  transport  including 
erosion/ entrainment,  diffusion,  settling,  deposition  and  consolidation.  Following  model 
testing  against  analytical  solutions,  laboratory  data  and  field  observations,  tidal  lutocline 
dynamics  was  examined  by  applications  of  COHYD-UF  and  COSED-UF  and  comparisons 
with  field  observations  on  flow  and  sediment  transport  from  the  Jiaojiang. 

7.2  Conclusions 

Important  conclusions  drawn  from  this  study  are  as  follows: 

1 . It  is  shown  that  turbulence  damping  in  the  water  column  at  high  SSC  is  governed 
by  the  settling  flux  in  the  flocculation  settling  range,  and  by  cohesion  and  interaction  between 
floes  in  the  hindered  settling  range.  Maximum  turbulence  damping  is  shown  to  occur  at  the 
lutocline,  supporting  a similar,  but  qualitative  observation  (e.g.,  Ross  and  Mehta,  1989).  Data 
derived  from  the  Jiaojiang  estuary  are  shown  to  support  this  observation. 

2.  Numerical  model  tests  using  SSC  data  from  the  Jiaojiang  showed  that  without 
considering  turbulence  damping,  the  simulated  results  of  SSC  become  less  tenable  when 
compared  with  the  field  observations,  resulting  in  comparatively  more  uniform  profiles  of 
SSC  and  less  distinct  lutoclines  than  observed. 
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3.  Observations  in  the  Jiaojiang  showed  that  the  lutocline  strength  was  highly 
associated  with  turbulence  energy.  The  lutocline  became  comparatively  less  distinct  as  the 
turbulent  energy  production  F >~0.0003  m ^ s This  supports  the  qualitative  conclusion 

from  the  derived  expression  of  turbulence  damping. 

4.  High  and  low  frequency  internal  waves  were  detected  at  the  lutocline  in  the 
Jiaojiang.  The  shallow  water  low  frequency  wave  had  a representative  rms  height  of  0.38  m 
and  a modal  frequency  of  0.09  rad  s * , which  was  near  the  local  Brunt-Vaisala  frequency. 
The  deep  water  high  frequency  wave  was  characterized  by  sharp  crests  and  flat  troughs,  with 
a rms  height  of  0.21  m and  a modal  frequency  of  1.33  rad  s '* . This  wave  was  possibly 
induced  by  interfacial  shear  at  the  lutocline. 

5.  The  height  and  the  angular  frequency  of  both  high  and  low  frequency  waves 
decreased  with  increasing  Richardson  number.  The  height  and  angular  frequency  versus 
logi?/Q  plots  exhibited  linear  trends  and,  in  turn,  the  celerity  and  wave  length  increased  with 

increasing  Richardson  number. 

6.  Numerical  modeling  tests  showed  that  the  developed  codes,  COHYD-UF  and 
COSED-UF,  are  able  to  adequately  simulate  transport  processes  including  sediment 
propagation,  advection,  diffusion,  entrainment  and  consolidation.  The  simulated  results 
compare  reasonably  with  laboratory  and  field  observations,  once  the  appropriate  process- 
related  parameters  are  adopted. 

7.  From  the  simulations  and  observations  in  the  Jiaojiang  estuary,  the  lutocline  was 
found  to  behave  as  follows:  (1)  The  lutocline  elevation  was  higher  during  peak  flows  than 
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during  slack  water,  with  the  highest  elevations  during  peak  floods;  (2)  the  lutocline  layer  was 
thinner  during  slack  water  than  during  peak  flows,  and  also  thinner  during  neap  tide  than 
spring;  (3)  the  lutocline  layer  rose  during  flow  acceleration  and  fell  during  flow  decelerating 
periods;  and  (4)  the  lutocline  was  observed  to  persist  through  most  of  the  tidal  cycle,  except 
1-2  hr  around  high  water.  The  overall  behavior  of  lutocline  reflects  the  cumulative  effects 
of  tidal  current  asymmetry  (with  a stronger  flood  peak  and  a weaker,  longer  duration  of  ebb 
peak),  sediment  settling  and  entrainment,  turbulence  damping  and  also  the  large  ratio  of  tidal 
range  to  mean  water  depth  resulting  in  a dilution  effect  during  flood. 

7.3  Recommendations  for  Future  Studies 

The  present  study  is  based  on  the  limited  data  from  the  Jiaojiang.  There  persists  a 
lack  of  quantitative  information  on  internal  wave  generation,  propagation,  interfacial 
instability  and  vertical  mixing  at  the  lutocline.  The  only  observations  used  to  verify  the 
derived  analytical  model  of  mixing  length  in  the  sediment-stratified  water  column  were  the 
results  obtained  from  the  vertical  excursion  of  SSC  using  ASSM  signals. 

In  the  sub-model  for  sediment  erosion,  it  was  assumed  that  the  existence  of  fluid  mud 
has  no  tangible  effect  on  bottom  erosion.  However,  in  general  fluid  mud  appears  as  a 
“protective”  cover  over  the  bottom  and  thereby  bottom  erosion. 

Based  on  the  conclusions  of  this  study,  the  following  recommendations  are  made  for 
future  research: 

1 . In  order  to  fully  understand  the  internal  wave  behavior,  collecting  long  time  series 
of  continuous  in  situ  records  of  the  interface  are  essential. 
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2.  Direct  measurement  of  turbulent  mixing  in  the  fluid  mud  is  necessary  to  verify  the 
mixing  length  model. 

3.  Extensive  experiments  in  flumes  are  required  for  understanding  the  way  in  which 
the  fluid  mud  layer  affects  bottom  erosion. 


APPENDIX  A 

DERIVATIONS  OF  THE  GOVERNING  EQUATIONS 

A.l  Vertical  Velocities,  w and  (o.  and  Continuity  Equation  Q.D 
The  o-transform  converts  the  time  and  Cartesian  coordinate  system  (t,  x,  y,  z)  to  the 
new  time  and  coordinate  system  (/',  x',y',  a)  according  to 

t'=t 


x'=x 


y'=y 


(A.1) 


H 


The  vertical  velocity  in  the  o-coordinate,  w,  is  defined  as 


(0= — 
Dt 


(A.2) 


The  vertical  velocity  in  the  z-coordinate,  w,  is 


w= — =H — +o 
Dt  Dt 


Do  DH  DC 
— +o + — - 


Dt  Dt  Dt 


(A.3) 


The  continuity  equation  is 


du  0v  dw 
— + — + — 


(A.4) 


dx  dy  dz 


The  transformation  of  each  term  in  Eq.  (A.4)  is  as  follows 


I8I 


182 


The  first  term: 


— da  _ du  _du'  a 8H ^ 1 5C ' 

dx  dx'  dx  da  dx  dx'  da  ^ H dx'  H dx\ 


(A.5) 


The  second  term: 


The  third  term: 


dv  _ dv  dy'  ^ dv  do  _ 5v  dv 
dy  dy'  dy  da  dy  dy'  da 


a dH  ^ 1 ac 


Hdy'  Hdy'^ 


(A.6) 


dw  _dw  da  _ 1 dw 
dz  da  dz  H da 


(A.7) 


Differentiating  Eq.  (A.3)  vsdth  respect  to  o,  one  obtains 


aw  jjdiji  du 
— -H — + — 

da  da  da 


dH  ac 

o — +— 
dx'  dx' 


dH  dv 
+u + 

dx'  da 


dH  aci  dH  dH 
a — + — — \ +v — + 


, dy'  dy' J dy'  dt' 


(A.8) 


Substituting  Eqs.  (A.5)-(A.8)  into  Eq.  ( A.4),  simplifying  and  eliminating  the  superscript 
(prime)  for  convenience,  one  obtains 


ac  duH  dvH  „ao)  „ 
—5.  + + +H — =0 

dt  dx  dy  da 


(A.9) 


Integrating  Eq.  (A.9)  with  o from  -1  to  0 and  considering  the  facts  that  -0,  the 

continuity  equation  (3.1)  becomes 


dHu  dHv],  „ 

+ k/o=0 

a^:  dy  ) 


(A.10) 


A.2  Momentum  Equations  (3.2)  and  OJ) 


The  derivation  of  the  momentum  Equation  (3.2)  in  the  x-direction  is  given  here.  The 
derivation  of  equation  (3.3)  in  the  >^-direction  is  analogous.  The  momentum  equation  in  the 
x-direction  is 


du  ^ du  du  du  r 1 dp 

— +u — +v — +w — -fv=—^  + 

dt  dx  dy  dz  p 
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Here  the  rheological  effect  (the  last  term  on  the  right  hand  side)  of  fluid  mud  is  simply 

considered  by  assuming  that  its  rheological  behavior  satisfies  the  Bingham  model  (Odd  and 
Cooper,  1989). 

The  o-transform  of  terms  on  the  left  hand  side  of  Eq.  (A.l  1)  is 
The  first  term: 


du  _ du  dt'  ^du  da  _ du  du 

[ o a//^  1 ac) 

dt  dt'  dt  da  dt  dt'  da 

i Hdt'  * Hdt'  ^ 

(A.12) 


The  third  term: 


du  _ du  dy'  ^du  da  _ du  du(  a dH  1 aC  ' 
dy  dy'  dy  da  dy  dy'  da[Hdy'  Hdy'j 


(A.  13) 


The  fourth  term: 


^_du  da  _ 1 du 

dz  da  dz  H da  (^-14) 


Assuming  pressure  to  be  hydrostatic,  the  pressure  term  in  right  side  of  Eq.  (A.l  1)  becomes 
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where  is  the  fluid  density  at  the  water  surface.  Thus  the  a-transform  of  the  pressure  term 
is 


1 dp_  89^ 
p dx  p 


acax'^ac  aoVgr 

oJ 


ax'  dx  do  dx' 


dx'  dal  Hdx'^Hdx', 


Hda 


p dx  p J dx'  pdx'Jda  pdx'J  da 
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The  transformation  of  the  turbulent  diffusion  term  in  the  x-direction  is 


d^u 

_ a 

du 

_ du( 

ax  2 

dx 

dx' 

da 

//dx'^Ndx'l 


=-^^+Higher  Order  Terms 
ax'2  dx'^ 


Similar  to  Eq.  (A.  17),  the  turbulent  diffusion  term  in  the  >^-direction  is 


-^=-^^+Higher  Order  Terms=-^^ 
ay  2 a/2  ^,2 


The  turbulent  diffusion  term  in  the  vertical  direction  has  the  form 


_a 

=1A 

du 

dz 

i ’azj 

Hda 

^ H da  J 

(A.  17) 


(A.  18) 


(A.  19) 


Substituting  Eqs.  (A.3),  (A.5)  and  (A.12)-(A.19)  into  Eq.  ( A.ll),  simplifying  and 
eliminating  the  superscript  (prime)  for  convenience,  the  momentum  equation  (3.2)  in  the  x- 
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direction  is  obtained  as 
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A.3  Sediment  Conservation  Equation  (3. 151 
The  sediment  conservation  equation  is 

dc  dc  dc  dc  ^ 


-+M +V +W- 


dt  dx  dy  dz  dz 


dx^  ^dy^ 


dz 


K — \ 
''dz 


Similar  to  Eq.  (A.20),  the  following  forms  are  obtained: 
The  first  term  on  the  left  hand  side  (LHS): 


dc  dc  dt'  dc  da  dc  dc 


dt  dt'  dt  da  dt  dt'  da 


a dH  I dC 


Hdt'  Hdt' 


The  second  term  on  the  LHS: 


dc  _ 
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The  forth  term  on  the  LHS 
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dc  _dc  do  _ 1 dc 
dz  do  dz  H do 


(A.25) 


The  fifth  term  on  the  LHS: 


3o)  c 1 5(0  c 

J _ 1 S 

dz  ~~H  do 


The  first  term  on  right  hand  side  (RHS): 

5^c  5^c 


dx~^  dx'^ 

The  second  term  on  the  RHS: 
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The  third  term  on  the  RHS: 
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(A.27) 


(A.28) 


(A.29) 

Substituting  Eqs.  (A.22)-(A.29)  into  Eq.  (A.21),  simplifying  and  eliminating  the  superscript 
(prime)  for  convenience,  the  sediment  conservation  equation  (3.15)  becomes 

dc  dc  dc  ( 

■+u — +v +(0- 
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APPENDIX  B 

NUMERICAL  TECHNIQUES 


B.I  Back-Tracing  Approach 

In  the  Eulerian-Lagrangian  differential  scheme,  the  physical  properties  of  water 
particles  at  time  step  n are  obtained  by  the  back-tracing  approach  (Casulli  and  Cheng,  1992). 
The  back-tracing  time  interval  At"=At/N^.  At  each  back-tracing  step  m,  the  water  particle 

will  lag  by  small  distances 


dx"'  = -u"  — 

AX  Ay 


da'"=-w‘ 


n At" 


AO 


(B.I) 


where  dx dy  ""and  do"’  are  the  backward  distances  at  time  step  m in  the  x,  and  o 
directions,  respectively,  and  u^,  and  (o^  are  the  local  velocities  of  the  water  particle  in 

the  X,  y and  a directions,  respectively,  that  are  approximated  by  bilinear  interpolation  over 
the  eight  surrounding  mesh  points.  Finally,  the  water  particle  is  traced  back  to  the  position 
p,  using  back-tracing  distances  dx,  dy  and  do  given  by  (Figure  B.I) 


dx  = ^ dx  dy='^  dy  do  = '^  do"" 

m = I 


m = l 


(B.2) 


m = l 


Then,  the  physical  property  at  point  p is  obtained  by  bilinear  interpolation  as  follows 
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Gp={\ -da){\  -dx){\ -dy)G^^^{\ -dx)dyG2 +dx{\ -dy)G"+dxdyG” 
+i/o[(l  -dx){\  -dy)Gs+(\  -dx)dyG^"+dx(\  -dy)G"+dxdyG” 


(B.3) 


Figure  B.l.  Schematic  diagram  of  back-tracing  approach,  where  dotted  line  is 
the  pathline  of  water  particle,  o is  the  position  of  water  particle  at  current  time 
step  n+1  and  p is  the  position  of  water  particle  at  the  previous  time  step  n. 


B.2  Pre-conditioned  Conjugate  Gradient  Method 
The  linear  five-diagonal  system  of  equations  for  the  water  surface  elevation, 
C,  can  be  written  in  the  following  general  form 


r - 


a. 


^ij*\  ^ij 


(B.4) 
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where  a.^j/2  and  b..  are  coefficients  that  are  dependent  on  the  time  step.  Note 

that  for  notational  simplicity,  superscript  («+l)  of  water  surface  elevation,  C,  has  been 
omitted.  Note  also  that  the  coefficients  and  cijj^y2  non-negative  and  their 

sum  is  strictly  less  than  unity.  Thus  the  system  formed  by  these  equations  is 
normalized,  symmetric  and  positive-definite  (Casulli  and  Cheng,  1992). 

The  pre-conditioned  conjugate  gradient  algorithm  to  solve  the  system  of 
equations  (B.4)  takes  the  following  steps  (Bertolazzi,  1990):. 


(B.2.1)  Guess  c!?. 


(B.2.2)  Set 


, -A 


(B.2.3)  Carry  out  following  calculations  for  k=0,  1,  2,  and  until 


(^(9,  ^(*)) 


(B.5) 
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where  e is  the  allowed  error.  At  each  iteration  the  essential  calculations  consist  of  a 
matrix-vector  multiplication  Mp^  ^ as  specified  in  Eq.  (B.5),  two  scalar  products 

between  vectors,  namely  and  and  three  sums  between  vectors, 

namely  and 


APPENDIX  C 

EFFECT  OF  TEMPERATURE  ON  SETTLING  VELOCITY 


The  experimental  results  of  Lau  (1994)  related  to  temperature  effect  on  the 
deposition  of  cohesive  sediments  in  an  annular  flume  indicated  that  temperature  affects  the 
settling  velocity  not  only  through  changes  in  the  viscosity  and  density  of  water,  but  also 
through  floe  aggregation.  The  floes  become  stronger  and  denser  with  decrease  of 
temperature,  and  consequently  their  settling  velocity  increases.  To  examine  this  effect,  Lau’s 
data  are  reprocessed  here.  His  experiments  were  carried  out  at  different  temperatures  and 
under  the  same  flow  condition  with  bottom  shear  stress  t^=0.2  Pa  and  initial  vertical  mean 

SSC,  Cq=9  kg  m . The  available  data  are  the  concentration-time  curves  for  runs  at  various 


temperatures.  Commercial  kaolinite  was  used  in  the  experiments.  Figure  C.l  shows  the 
frequency  distribution,  4y,  of  the  settling  velocity  based  on  the  standard  size  distribution  of 


a similar  kaolinite  (Yeh,  1979),  where  the  settling  velocity  is  calculated  from  the  Stokes  law 


0« 


2( 


18v 


P. 


— -1 

iPo 


(C.l) 


In  Eq.  (C.l),  and  o)q^,  respectively,  are  the  grain  size  and  settling  velocity  of  the 
kaolinite  class,  and  =2,650  kg  m^  po=l,000  kg  m^and  v = 10‘®  m^  s"‘. 
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In  accordance  with  the  method  of  Mehta  and  Lott  (1987),  sorting  by  size 
characteristically  occurs  for  non-uniform  fine  sediment.  Given  C(t)  as  the  instantaneous 
depth-averaged  concentration,  under  a given  depositional  flow  condition  C will  decrease 
with  time  and  ultimately  reach  a steady  state  value,  Cj.  If  and  are  respectively 

defined  as  the  minimum  and  maximum  critical  shear  stress  for  deposition  of  the  m ^ kaolinite 

classes,  then  when  no  sediment  can  deposit,  whereas  when  , all  sediment 

eventually  deposits.  When  a certain  fi-action  of  the  initial  sediment,  represented 

by  C^,  will  ultimately  remain  in  suspension,  and  the  remainder,  represented  byCg-C^  will 

settle  out.  The  occurrence  of  Cj  less  than  Cq  but  greater  than  zero  is  an  indication  of 


sediment  sorting.  A deposition  law  for  a non-uniform  cohesive  sediment  was  developed  by 
Mehta  and  Lott  (1987)  based  on  the  consideration  that  the  instantaneous  concentration,  C, 
is  obtained  by  summation  of  the  corresponding  concentrations,  C^,  obtained  fi-om  the 


deposition  relationship  of  Krone  (1962)  for  each  class.  This  leads  to 


C(0  _ 1 


M, 


M. 


E C„(0=E  4>/wJexp^ 


^0  ^0 

where  is  the  settling  velocity  of  the  kaolinite  class,  which  is  simplified  as 


1— ^ 

( \ 

0) 

sn  ^ 

H 

(C.2) 


(C.3) 

where  and  are  respectively  the  minimum  and  maximum  settling  velocities  of  the 

kaolinite  classes  and  is  a flocculation  factor. 
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Using  Eqs.  (C.2)  and  (C.3),  one  can  attempt  to  best-fit  the  experimental  data  of  Lau 
(1994)  for  different  temperatures  by  changing  the  minimum  settling  velocity,  , and  the 

flocculation  factor,  F^.  The  results  are  shown  in  Figures  (C.2)-(C.5),  where 

taken.  The  plot  of  cumulative  weight  finer  against  the 

settling  velocity  for  different  temperatures  is  shown  in  Figure  (C.6).  It  is  observed  that  the 
settling  velocity  increases  with  decreasing  of  temperature.  In  other  words,  the  floes  become 
stronger  and  denser  with  decreasing  of  temperature. 

To  express  the  temperature  effect  mathematically,  the  (median)  settling  velocity  at 
50%  of  cumulative  weight  (finer)  is  taken.  The  results  are  shown  in  Figure  (2.3),  where  the 
dimensionless  median  settling  velocity  is  plotted  against  the  temperature  with  as  the 

settling  velocity  at  15  °C.  Finally,  from  Figure  (2.3)  the  temperature  function,  F,  i.e.,  Eq. 
(2.28),  is  obtained  as 

1.776 -0.05 186,  for  0=0-30  °C 


(C.4) 


Frequency  distribution, 
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Figure  C.l.  Frequency  distribution,  (|y,  of  the  settling  velocity  of  kaolinite. 


Figure  C.l.  Time-concentration  relationship  during  deposition  at 
26  °C.  Open  circles  are  the  experimental  data  of  Lau  (1994). 


Normalized  concentration  C/C„  Normalized  concentration  C/C 
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Time  Chr) 

Figure  C.3.  Time-concentration  relationship  during  deposition  at 
20  °C.  Open  circles  are  the  experimental  data  of  Lau  (1994). 


Figure  C.4.  Time-concentration  relationship  during  deposition  at 
10  °C.  Open  circles  are  the  experimental  data  of  Lau  (1994). 
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Figure  C.5.  Time-concentration  relationship  during  deposition  at 
5 °C.  Open  circles  are  the  experimental  data  of  Lau  (1994). 


Figure  C.6.  Cumulative  distribution  of  settling 
velocity  of  kaolinite  at  different  temperatures,  0. 


APPENDIX  D 

AN  APPLICATION  OF  COHYD-UF:  CONTRACTION  SCOUR  IN  A RIVER 

D.l  Scour  Problem 

An  application  of  COHYD-UF  is  reported  here  for  the  simulation  of  a contraction 
scour  problem  at  the  Haldia  Jetty,  or  Pier  (Figure  D.  1 ),  in  the  Hooghly  River  in  India  (Figure 
D.2).  The  pier,  designed  in  the  1960s,  served  as  an  oil  unloading  terminal  that  enabled 
40,000  DWT  oil  tankers  to  transfer  oil  to  storage  tanks  at  the  Haldia  Port.  The  Hooghly 
estuary  transports  high  loads  of  fine-grained  material.  Since  constmction  of  the  pier,  a scour 
occurred  in  front  of  the  pier.  The  presence  of  the  hole  near  the  tip  of  the  pier  was  detrimental 
to  the  pier  piles  and  hence  the  superstructure  (Engineers  India  Limited,  1980;  Rao,  et  al., 
1980).  The  surficial  shape  of  the  scour  hole  was  nearly  elliptical,  with  the  major  axis  and 
minor  axis  approximately  perpendicular  and  parallel,  respectively,  to  the  longitudinal  axis 
of  the  pier.  The  scour  depth  varied  between  7 m and  13  m,  and  the  hole  migrated  southward 
over  a 40-month  period  up  to  May,  1980.  The  contours  in  Figure  D.3  show  the  hole  as  it 
appeared  in  the  May,  1980  survey.  Note  that  the  scour  depths  are  taken  as  those  below  the 
original  bottom  elevation  before  the  construction  of  the  pier. 

During  the  40-month  period  the  average  rate  of  migration  was  1 .53  m/month  along 
the  direction  of  the  main  flow  and  0.56  m/ month  perpendicular  to  the  main  flow  and  away 
from  the  pier.  This  mode  and  rate  of  migration  suggested  a slight  dominance  of  the  ebb 
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current  over  flood  in  the  bottom  region  immediately  adjacent  to  the  pier.  The  strength  of  ebb 
current  in  this  region  was  on  the  order  of  2 m s ‘.  The  presence  of  such  a strong  current, 
coupled  with  the  fact  that  a 40,000  DWT  tanker  was  often  docked  at  the  pier,  exacerbated 
the  scour  problem  and  necessitated  action  for  filling  up  the  hole  in  order  to  stabilize  the  pier 
(Engineers  India  Limited,  1980). 


Figure  D.  1 . Schematic  diagram  showing  the  Haldia  oil  pier  and  depth 
contours  (m)  in  the  vieinity.  Water  depth  are  below  mean  low  water. 

An  evident  conclusion  that  can  be  drawn  is  that  the  reduetion  of  the  flow  area  due  to 
the  pier  resulted  in  an  increase  in  the  high  ebb  current  and  associated  bed  shear  stress  in  front 
of  the  pier.  Consequently,  there  was  an  inerease  in  the  erosive  force  in  the  effectively 
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contracted  area  in  front  of  the  pier,  which  is  believed  to  be  the  main  cause  of  the  scour  hole. 


Hence,  it  can  be  called  contraction  scour  (HEC-18,  1995). 


Calcutta 
Garden  Reach 


Study  area 


Beaumont  % 
Gut 


Figure  D.2.  Location  map  of  Haldia  oil  pier,  India. 
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D.2  Scour  Simulation 

As  described  above,  the  disturbance  of  the  pier  on  the  local  flow  and  the  dominance 
of  ebb  current  were  thought  to  be  the  main  reasons  for  the  development  of  the  hole. 
Accordingly,  the  scour  depth  can  be  calculated  through  the  modeled  ebb  flow,  which  will  be 
assumed  to  be  constant.  Figure  D.4  shows  the  bathymetry  within  the  modeled  segment  of  the 
river.  The  bed  erosion  rate  is  taken  to  be  proportional  to  the  bottom  shear  stress  t.  as 

follows 


m=M  ^ 

e max 


(D.l) 


\ ^ y 

By  assuming  the  bottom  to  be  in  sedimentary  equilibrium  before  the  construction  of  the  pier, 
one  can  take  the  critical  shear  stress  for  erosion,  to  be  the  bottom  shear  stress  before  pier 


construction.  The  erosion  rate  constant,  can  be  estimated  using  the  maximum  scour 


depth,  a//^,  according  to 


aH 
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(D.2) 


where  x^'  and  x^'  respectively  are  the  bottom  and  critical  shear  stresses  at  the  point  where 
the  maximum  scour  occurred,  aT  is  the  time  period  over  which  a//^^  occurred,  and  is 
the  dry  density  of  the  bottom  sediment.  Accordingly,  the  scour  depth  at  any  position,  aH,  is 


obtained  from 
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Figure  D.3.  Measured  scour  depths  in  front  of  the  Haldia  oil  pier. 
The  pier  is  shown  as  an  idealized  rectangular  protrusion.  Unit:  m. 


East  Bank  of  River 


Figure  D.4.  Bottom  topography  of  the  modeled  segment  of  the  river  in  the 
vicinity  of  the  Haldia  pier.  Water  depths  (unit:  m)  are  below  mean  low  water. 
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a//= 


\ 


(D.3) 


where  the  bottom  shear  stress,  was  calculated  from  Eq.  (2.9)  without  considering 
stratification  effects.  The  values  of  Zq=0.1  mm  [effective  roughness  of  the  bed  in  Eq.  (2.9)] 
^ taken  in  these  calculations.  From  the  modeling  tests,  it  was  found  that 
the  ratio  of  t,'  to  x 'was  about  1.51. 

D S 

The  modeled  region  includes  the  area  from  2.25  km  upstream  to  2.25  km  downstream 
of  the  pier;  thus  with  a length  of  4.5  km  (Figure  D.4).  The  domain  was  discretised  with 
spatial  steps  ax=50  m,  a>'=25  m and  ao=0.1,  with  the  total  number  of  rectangular  cells 
MxiVxZ^=90x  104x10  =93,600.  A time  step  a/=5  s was  adopted  due  to  stability  constraints 

resulting  from  the  numerical  scheme  involving  in  the  back-tracing  approach  mentioned  in 
Chapter  2 [Eq.  (2.50)].The  lowest  layer  near  the  bottom  was  0.05// above  the  bed,  and  the 
highest  layer  near  the  surface  was  0.05// below  the  local  water  surface. 

The  upstream  and  downstream  open  boundary  conditions  were  prescribed  by  constant 
water  surface  elevations  at  both  ends.  Elevations  of  0.25  m above  mean  low  water  level  at 
the  upstream  boundary  and  0.00  m at  downstream  open  boundary  were  taken.  Model  run  was 
initiated  at  zero  velocity  and  a constant  water  surface  slope  of  5.56x10'^  over  the  domain. 


COHYD-UF  was  run  until  a stable  ebb  flow  field  resulted.  The  model  was  run  for  two  cases: 
before  and  after  the  construction  of  the  pier,  and  both  results  were  outputted  for  above 
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calculations.  At  each  grid  point,  the  bottom  shear  stress,  in  Eq.  (D.3),  was  taken  as  the 
simulated  value  after  the  construction  of  the  pier,  and  the  bottom  shear  strength,  x^  in  Eq. 
(D.3),  was  taken  as  the  x^  value  before  the  construction  of  the  pier. 

D.3  Results 

Figures  D.5  and  D.6  show  simulated  flow  fields  around  the  pier  near  the  surface  and 
the  bottom,  respectively.  It  is  observed  that  the  current  in  the  front  of  the  pier  became 
stronger  due  to  the  disturbance  of  the  pier  on  flow,  with  a maximum  ratio  between  velocities 
after  and  before  the  construction  of  the  pier  of  about  1 .23. 

The  solid  contours  in  Figure  D.7  show  the  simulated  scour  depths  in  the  front  of  the 
pier.  In  Figure  D.8,  the  simulated  areas,  A^,  corresponding  to  2-3  m,  3-4  m,  4-5  m and  >5 

m scour  depths  are  plotted  against  the  corresponding  measured  areas,  A^ . This  comparison 

shows  that  the  model  reasonably  reproduced  the  area  of  the  scour  hole,  and  that  the  deeper 
the  scour  depth,  the  better  the  comparison  of  the  simulated  versus  measured  scour  area. 


Distance  (m)  Distance  (m) 
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Figure  D.5.  Simulated  flow  field  around  pier  at  0.05// below  the  surface. 


Distance  (m) 


Figure  D.6.  Simulated  flow  field  around  pier  at  0.05// above  the  bottom. 


Simulated  scour  area,  (m^) 
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Figure  D.7.  Comparison  of  scour  depths  simulated  (solid  lines) 
and  measured  (dashed  lines)  in  front  of  the  Haldia  oil  pier.  Unit:  m. 


Figure  D.8.  Comparison  between  simulated  and  measured  areas 
at  2-3  m (•),  3-4  m (^),  4-5  m (O)  and  >5  m (+)  scour  depths. 


APPENDIX  E 

SIMULATION  OF  SEDIMENT  DEPOSITION  IN  A FLUME 

E.l  Introduction 

In  order  to  test  the  ability  of  COSED-UF  in  predicting  sedimentation,  a simulation 
was  conducted  to  compare  modeled  and  observed  shoaling  patterns  in  the  single  flume  test 
of  Ariathurai  (1974).  In  the  simulation,  the  settling  velocity  in  moving  water  was  determined 
using  experimental  data  of  cohesive  sediment  deposition  from  Mehta  (1973). 

E.2.  Flume  Test 

In  the  test  of  Anathurai  (1974),  a 20  m long,  61  cm  wide  tilting  recireulating  flume 
was  used.  Approximately  two-thirds  of  the  way  down  the  flume,  a grating  made  of  vertieal 
steel  rods,  each  of  6.35  mm  diameter,  was  plaeed  across  one-half  of  the  flume  width  so  that 
it  partially  blocked  the  flow.  The  barrier  allowed  eonsiderable  flow  through  itself 
Reeonstituted  seawater  with  a salinity  of  about  3 1 %o  was  added  to  the  flume  to  give  a flow 
depth  of  10  cm.  The  flume  pump  was  then  set  so  as  to  generate  an  average  veloeity  of  about 
1 7 cm  s . The  slope  of  the  flume  bottom  was  adjusted  simultaneously  to  produee  a uniform 
flow  depth  along  the  length. 

Before  starting  the  experiment,  the  barrier  was  removed  and  the  flow  velocity 
increased.  Sediment  was  then  added  to  the  flume  gradually.  San  Francisco  Bay  sediment  (bay 
mud)  obtained  from  deposits  in  a yaeht  harbor  in  Mare  Island  was  used  for  the  test.  Sea 
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shells,  silt  and  other  coarse  materials  present  in  the  sample  were  removed  by  sedimentation 
and  the  remaining  cohesive  sediment  with  some  silt  in  it  was  mixed  in  sea  water  and  poured 
into  running  water  in  the  flume  slowly  until  an  initial  mean  SSC  of  0.2  kg  m was  reached. 

After  all  of  the  sediment  had  been  added  the  barrier  was  placed  and  the  velocity 
decreased  to  the  same  value  as  earlier  (i.e.,  17  cm  s ’' ).  Afterwards,  the  suspended  sediment 

settled  gradually  in  the  region  in  the  lee  of  the  barrier  having  relatively  low  velocities.  The 
test  was  earned  out  for  a period  of  about  3.5  days.  Water  was  then  drained  out  slowly  and 
the  deposition  pattern  around  the  barrier  was  recorded. 

It  was  found  that  the  SSC  decreased  with  time  according  to  (Ariathurai,  1974) 

C-CqIO  (E.l) 

where  C is  the  vertically-averaged  SSC,  Cg  is  the  initial  SSC  (==0.2  kg  m and  is  a 

deposition  rate  constant,  which  was  found  to  be  1.0x10"^  s by  calibration  (Krone,  1962; 
Ariathurai,  1974). 

E.3  Settling  Velocity  in  Moving  Water 

As  stated  in  Section  4.3.5,  the  settling  velocity  in  moving  water  is  characteristically 
different  from  that  in  quiescent  water,  since  increasing  turbulence  can  enhance  flocculation 
and  at  the  same  time  limit  the  size  of  floes  that  can  be  sustained.  To  determine  the  relevant 
settling  velocity,  experimental  data  on  bay  mud  deposition  reported  in  Mehta  (1973)  were 
used.  These  experiments  were  carried  out  in  flumes  under  different  bottom  shear  stress,  t,. 
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and  initial  SSC,  Cq  . The  data  have  been  reported  as  the  fraction  of  depositable  concentration, 
C * [=(Cq-C)/ (Cq-C^)]  , as  a function  of  non-dimensional  time,  t * where  Cj  is 

the  final  (steady  state)  SSC  and  is  the  time  corresponding  to  C ’=50%.  Table  E.l  gives 
the  basic  parameters  in  these  experiments. 


Table^BJ_^_Basic^p^ameter^^  experiments  using  the  bay  mud 


No. 

Co 

(kg  m'3) 

H 

(m) 

(Pa) 

(Pa) 

9 

(kg  m-5) 

^50 

(hr) 

Investigator 

1 

8.45 

0.152 

0.165 

0.194 

4.23 

10 

Mehta  (1973) 

2 

0.50 

0.305 

0.067 

0.070 

0.05 

63 

Krone (1962) 

3 

0.92 

0.305 

0.049 

0.065 

0.00 

23 

Krone  (1962) 

4 

1.92 

0.305 

0.037 

0.065 

0.00 

7 

Partheniades 

(1962) 

5 

21.00 

0.152 

0.031 

0.065 

0.00 

5.9 

Krone  ( 1 962) 

Once  C {t  ),  /jQ,  Cq,  Cj  and  water  depth  H are  known,  the  settling  velocity  can  be 


calculated  from 


) 


Ca(p, 


(E.2) 


where  C”  and  C"’*  are  SSC  at  two  consecutive  measurements,  a/  is  the  time  interval 
between  these  two  measurements,  is  the  probability  of  sediment  deposition,  which  is 


209 


taken  as  1 x^/x^  (Krone,  1 962),  x^  is  the  critical  shear  stress  for  deposition,  and  the 

instantaneous  concentration  C=(C  ”+C"^‘)/2 . When  applying  Eq.  (E.2),  one  should  consider 

the  non-uniformity  of  the  sediment.  In  this  case  each  class  of  the  sediment  has  a different 
value  of  (Appendix  C).  By  assuming  that  there  is  no  interaction  among  different  sediment 

classes,  Mehta  and  Lott  (1987)  related  the  critical  stress,  x^^,  to  the  settling  velocity, 

for  the  n sediment  class  by 


where  [ -ln(T^J^^/T^,)/ln(a)^J^/(o^,)]  is  a sediment-dependent  constant,  e.g.,  p^=0.5  for 

kaolinite  (Mehta  and  Lott,  1987).  The  other  symbols  are  the  same  as  in  Appendix  C.  For  the 
bay  mud,  t^j=^0.065  Pa,  t^j^^=0.318  Pa,  and  2.3x10'^  ms"*  were  chosen  (Mehta, 

1973).  Considering  that  under  a given  flow  only  sediment  classes  having  x,  are 

® dn  b 

depositable,  one  can  take  x^  in  Eq.  (E.2)  as  the  minimum  value  of  the  critical  stress  of  the 
depositable  sediment  classes,  i.e.,  t^=MIN(t^Jt^^^Tj).  Table  E.l  gives  for  each  test. 

The  calculated  results  using  Eq.  (E.2)  are  shown  in  Figure  E.l,  where  the  solid  line 
is  the  best- fit  of  the  calculated  data  points  using  the  settling  velocity  - SSC  relation  given  by 
Eq.  (4. 12).  It  is  seen  that  for  SSC<~0.2  kg  m the  particles  become  practically  free  settling, 
with  a settling  velocity  of  about  4.0  x 10"^  m s . 
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Figure  E.l.  Settling  velocity  as  a function  of  SSC 
in  moving  water.  Data  are  from  Mehta  (1973). 


E.4  Deposition  Simulation 

In  the  simulation,  the  flume  was  discretized  with  spatial  steps  ax=0.26  m, 
Ay=0.02032  m and  ao=0.1,  with  the  total  number  of  rectangular  cells  A/XiVxZ^=77x30xio 

=23,100.  A time  step  a/=0.02  s was  adopted  due  to  the  same  reason  as  stated  in  Appendix 
D.  The  lowest  layer  near  the  bottom  was  0.057/ above  the  bed  and  the  highest  layer  near  the 
surface  was  0.05// below  the  local  water  surface. 
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The  upstream  and  downstream  boundary  conditions  were  prescribed  by  constant 
water  surface  elevations  at  both  ends.  Elevations  of  2.5  cm  below  the  horizontal  datum  at  the 
downstream  and  0.00  m at  the  upstream  boundary  were  taken.  This  yielded  a mean  water 
surface  slope  of  1.25x10  ^ . The  flume  bottom  was  tilted  at  the  same  slope  as  the  water 
surface.  The  water  depth  was  10  cm.  The  following  values  of  the  hydrodynamic  parameters 
were  empirically  selected:  ^^=0.05  m^  s"’,  ^^=0.005  m^  s'*  and  z^=0.5  mm.  Under 

these  conditions,  the  model  generated  a cross-sectional  mean  velocity  of  17  cm  s before 
the  setting  of  the  barrier. 

In  the  model  grid,  the  real  barrier  was  approximated  by  alternatively  placed  solid 
lines  at  the  same  location  as  in  the  flume  test.  The  model  was  initiated  at  zero  velocity  and 
a constant  water  surface  slope  of  1.25x10-3.  COHYD-UF  was  then  run  until  a stable  flow 

field  resulted.  The  stable  result  was  outputted  for  calculation  of  deposition.  Figure  E.2  shows 
the  simulated  flow  field.  It  is  observed  that  the  velocity  is  relatively  weak  in  the  regions 
before  and  below  the  barrier,  and  also  near  the  side  walls,  which  are  potential  areas  of 
deposition  (Ariathurai,  1974). 

As  described  earlier,  during  the  deposition  process  the  SSC  in  the  flume  was  found 
to  be  a function  of  time,  as  defined  by  Eq.  (E.l).  Thus,  the  deposition  depth  at  each  grid  point 
can  be  simply  evaluated  from 


^ Co)  At 


1— ^ 


(E.4) 
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where  is  the  dry  density  of  the  newly  deposited  sediment  (140  kg  m'^).  The 


instantaneous  concentration,  C,  was  calculated  by  Eq.  (E.l),  and  the  bottom  shear,  from 

’ o’ 


Eq.  (3.9). 


0 0.5  1 1.5  2 2.5  3 3.5 


(a) 


0 0.5  1 1.5  2 2.5  3 3.5 


Distance  (m) 


Figure  E.2.  Simulated  flow  field  around  the  barrier  in 
the  flume,  (a)  near  the  surface  and  (b)  near  the  bottom. 


The  calculated  sediment  deposition  in  the  region  downstream  the  barrier  are  shown 
in  Figure  E.3.  Also  shown  in  this  figure  are  the  measured  thickness  of  the  deposit.  It  is  seen 
that  the  majority  of  deposition  took  place  downstream  the  barrier  due  to  a dramatic 
decreasing  in  the  velocity  there  (Figure  E.2). 


Distance  (m) 
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Figure  E.3.  Distribution  of  simulated  (solid  lines)  and  observed  (numbers  in  circles) 
deposition  (thickness)  at  the  down  side  of  the  barrier.  Data  are  from  Ariathurai  (1974). 
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